Why is the solution a scalar function, when it's defined on a VectorFunctionSpace?

1

I have defined a VectorFunctionSpace V on a domain and u is the unknown displacement vector defined on V, which should be found out by minimizing a strain-energy function. My problem is in 3-D. So for a point on the domain (x_{1},y_{1},z_{1}), the displacement vector will be (u_{1},u_{2},u_{3}) = u(x_{1},y_{1},z_{1}). But on solving the problem, it turns out that u is a scalar.

For example consider the following FEniCS script followed by the gmsh script used to construct the geometry:

FEniCS script:

from dolfin import *     
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
import numpy as np
import matplotlib.pyplot as plt
from ufl import cofac

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh
    
#--------------------------------------------------------------IMPORTING GEOMETRY FROM GMSH----------------------------------------------------------------------------------------------
msh = meshio.read("MWE.msh")

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_custom = Measure("dS", domain=mesh, subdomain_data=mf) 

dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

#*****************************************************COMPILING OUTER 3 LAYERS TO SOLVE THE ELASTICITY PROBLEM***************************************************************************

combined_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
combined_subdomains.array()[cf.array()==5] = 1  # Assigning same number for all the submeshes for compiling them into one submesh
combined_subdomains.array()[cf.array()==6] = 1
combined_subdomains.array()[cf.array()==7] = 1
mesh_ima = SubMesh(mesh, combined_subdomains, 1)

boundary_marker_ima = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim()-1, 0)
ncells0 = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim())
surface_marker_ima = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim(), 0)
vmap0 = mesh_ima.data().array("parent_vertex_indices", 0)
cmap0 = mesh_ima.data().array("parent_cell_indices", mesh_ima.topology().dim())

        
n = 0
for c in cells(mesh_ima):
  parent_cell = Cell(mesh, cmap0[c.index()])
  surface_marker_ima.array()[c.index()] = cf.array()[parent_cell.index()]
  for f in facets(parent_cell):
    for g in facets(c):
      g_vertices = vmap0[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        boundary_marker_ima.array()[g.index()] = mf.array()[f.index()]
    n=n+1

ds_ima = Measure("ds", domain=mesh_ima, subdomain_data=boundary_marker_ima) 
dx_ima= Measure("dx", domain=mesh_ima, subdomain_data=surface_marker_ima)


#------------------------------------------------------SOLVING HYPERELASTICITY PROBLEM ON THE OUTER 3 LAYERS-----------------------------------------------------------------------------

V = VectorFunctionSpace(mesh_ima, "Lagrange", 1) # Introducing a vector function space on intima, media and adventitia combined, i.e., ima submesh

n = FacetNormal(mesh_ima) # Defining the normal to the both the boundaries, when supplied with the bounday measure in the variational form, it automatically takes the normal on one of                              the boundaries according to the boundary measure supplied to the variational form

#------------------------------------------------------DEFINING FUNCTIONS NEEDED FOR THE VARIATIONAL FORM--------------------------------------------------------------------------------

du = TrialFunction(V)            # Incremental displacement
u  = Function(V)                 # Deformation vector
v  = TestFunction(V)             # Test function 

#-------------------------------------------------------------SETTING THE KINEMATICS-----------------------------------------------------------------------------------------------------
d = mesh_ima.geometry().dim()
I = Identity(d)             
F = I + grad(u)             
g_int = Identity(d)  
g_med = Identity(d)  
g_ad = Identity(d)  
C_int = (F*inv(g_int)).T*(F*inv(g_int))                   
C_med = (F*inv(g_med)).T*(F*inv(g_med))                   
C_ad = (F*inv(g_ad)).T*(F*inv(g_ad))                   
M = cofac(F)
P = Constant(5.0)
  

Ic_int = tr(C_int)     
Ic_med = tr(C_med)     
Ic_ad = tr(C_ad)     
J_int  = det(F*inv(g_int))
J_med = det(F*inv(g_med))
J_ad = det(F*inv(g_ad))
J_e = det(F)
C_lumen = F.T*F
Ic_lumen = tr(C_lumen)

# Elasticity parameters 
nu = 0.497   
mu1, mu2, mu3 = 27.9, 1.27, 7.56 
eta_int, eta_med, eta_ad = 263.66, 21.6, 38.57
beta_int, beta_med, beta_ad = 170.88, 8.21, 85.03
rho_int, rho_med, rho_ad = 0.51, 0.25, 0.55

# Stored strain energy density 
psi1 = (mu1/2)*(Ic_int - 3) + ((nu)*(mu1)/(1-2*nu))*(J_int-1)**2 - mu1*(ln(J_int)) + (eta_int/beta_int)*(exp(beta_int*((1-rho_int)*(Ic_int - 3)**2))-1)
psi2 = (mu2/2)*(Ic_med - 3) + ((nu)*(mu2)/(1-2*nu))*(J_med-1)**2 - mu2*(ln(J_med)) + (eta_med/beta_med)*(exp(beta_med*((1-rho_med)*(Ic_med - 3)**2))-1)
psi3 = (mu3/2)*(Ic_ad - 3) + ((nu)*(mu3)/(1-2*nu))*(J_ad-1)**2 - mu3*(ln(J_ad)) + (eta_ad/beta_ad)*(exp(beta_ad*((1-rho_ad)*(Ic_ad - 3)**2))-1)

  
Pi = (det(g_int))*(psi1*dx_ima(5)) + (det(g_med))*(psi2*dx_ima(6)) + (det(g_ad))*(psi3*dx_ima(7)) - (0.5)*dot(M*P*(-n), u)*ds_ima(1)

# Compute first variation of Pi (directional derivative about u in the direction of v)
F1 = derivative(Pi, u, v)

# Compute Jacobian of F
Ju = derivative(F1, u, du)

bcs = []

# Solve variational problem
solve(F1 == 0, u, bcs,
  form_compiler_parameters=ffc_options)
  
print(u.vector()[:], u.vector()[:].shape)

The gmsh script (named MWE.geo) is:

SetFactory("OpenCASCADE");

x1 = 0;
x2 = 2*Pi;
n =200;
lc1 = 0.08;
lc2 = 0.08;
lc3 = 0.08;
lc4 = 0.08;

//First interface Points and Lines for the endothelium
For i In {0:n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 1.7217+0.2858*Cos (x)-0.2452*Sin (x)+0.0204*Cos (2*x)+0.0345*Sin (2*x)-0.1876*Cos (3*x)+0.0348*Sin (3*x)-0.0358*Cos (4*x)+0.0329*Sin (4*x)-0.0101*Cos (5*x)+0.1119*Sin (5*x);

  pList1[i] = newp;

  Point (pList1[i]) = {r*Cos (x),r*Sin (x), 0}; //if you want the inner boundary mesh to be more refined change lc1 to a smaller number

EndFor

//Spline creates a continuous boundary using the points above
s1 = newreg;
Spline(s1) = {pList1[{0:n-1}],1};

//Second interface Points and Lines for the IEL
For i In {n:2*n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 2.4656+0.2869*Cos(x)-0.2907*Sin(x)-0.0445*Cos(2*x)-0.0129*Sin(2*x)-0.1528*Cos(3*x)+0.0131*Sin(3*x)-0.0457*Cos(4*x)-0.0665*Sin(4*x)+0.0186*Cos(5*x)-0.0081*Sin(5*x);

  pList2[i] = newp;

  Point (pList2[i]) = {r*Cos (x),r*Sin (x), 0};

EndFor

//Spline creates a continuous boundary using the points above
s2 = newreg;
Spline(s2) = {pList2[{n:2*n-1}],n+1};

//Second interface Points and Lines for the MEDIA-ADVENTITIA interface
For i In {2*n:3*n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 4.9314+0.1441*Cos(x)-0.2625*Sin(x)-0.1752*Cos(2*x)-0.2901*Sin(2*x)+0.1908*Cos(3*x)+0.0297*Sin(3*x)+0.2856*Cos(4*x)-0.2104*Sin(4*x)-0.0537*Cos(5*x)-0.1570*Sin(5*x);

  pList3[i] = newp;

  Point (pList3[i]) = {r*Cos (x),r*Sin (x), 0};

EndFor

//Spline creates a continuous boundary using the points above
s3 = newreg;
Spline(s3) = {pList3[{2*n:3*n-1}],2*n+1};

//Second interface Points and Lines for the ADVENTITIA outermost boundary
For i In {3*n:4*n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 7.9536+0.0443*Cos(x)-0.5983*Sin(x)-0.0620*Cos(2*x)+0.1528*Sin(2*x)+0.4119*Cos(3*x)-0.1249*Sin(3*x)-0.1130*Cos(4*x)-0.2179*Sin(4*x)+0.0910*Cos(5*x)+0.1176*Sin(5*x);

  pList4[i] = newp;

  Point (pList4[i]) = {r*Cos (x),r*Sin (x), 0};

EndFor

//Spline creates a continuous boundary using the points above
s4 = newreg;
Spline(s4) = {pList4[{3*n:4*n-1}],3*n+1};

//+
Physical Curve("C1", 1) = {1};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {3};
//+
Physical Curve("C4", 4) = {4};
//+
Curve Loop(1) = {2};
//+
Curve Loop(2) = {1};
//+
Plane Surface(1) = {1, 2};
//+
Curve Loop(3) = {3};
//+
Curve Loop(4) = {2};
//+
Plane Surface(2) = {3, 4};
//+
Curve Loop(5) = {4};
//+
Curve Loop(6) = {3};
//+
Plane Surface(3) = {5, 6};
//+
Physical Surface("S5", 5) = {1};
//+
Physical Surface("S6", 6) = {2};
//+
Physical Surface("S7", 7) = {3};
//+
Curve Loop(7) = {1};
//+
Plane Surface(4) = {7};
//+
Physical Surface("S8", 8) = {4};

The geometry along with the mesh for the above gmsh script looks like:

Screenshot (116)
Screenshot (116)1352×761 67 KB

I want the solution to be a vector function of the coordinates such that (u_{1},u_{2},u_{3}) = u(x_{1},y_{1},z_{1}). What am I doing wrong?

Any suggestions or feedback would be extremely helpful.

2

When you solve a Finite element problem, it being on a vector, tensor or mixed space, the data is always flattened. This is because you are solving a system on linear equations on the form Au=b,u is the vector containing the degrees of freedom.

To print it as a set of vectors, you can use the strategy shown in: How to get the solution obtained by fenics as a vector? - #7 by dokken

1 Like
3

Hi @dokken ,

Although the code from the site you suggested works for the MWE I supplied. But the same code produces the following error when I am restricting the solution to a particular boundary.

vector1[i] = u_e.sub(i, deepcopy=True).vector().get_local()
AttributeError: 'numpy.ndarray' object has no attribute 'sub'

Here is the extension of the problem in the FEniCS script I provided at the start. I am working with the same gmsh script I provided initially.

from dolfin import *     
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
import numpy as np
import matplotlib.pyplot as plt
from ufl import cofac

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh
    
#--------------------------------------------------------------IMPORTING GEOMETRY FROM GMSH----------------------------------------------------------------------------------------------
msh = meshio.read("matlabfenicslumintmedadcoarse.msh")

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_custom = Measure("dS", domain=mesh, subdomain_data=mf) 

dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

V_orig = VectorFunctionSpace(mesh, "P", 1)  # Defining VectorFunctionSpace on the whole mesh

#*****************************************************COMPILING OUTER 3 LAYERS TO SOLVE THE ELASTICITY PROBLEM***************************************************************************

combined_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
combined_subdomains.array()[cf.array()==5] = 1  # Assigning same number for all the submeshes for compiling them into one submesh
combined_subdomains.array()[cf.array()==6] = 1
combined_subdomains.array()[cf.array()==7] = 1
mesh_ima = SubMesh(mesh, combined_subdomains, 1)

boundary_marker_ima = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim()-1, 0)
ncells0 = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim())
surface_marker_ima = MeshFunction("size_t", mesh_ima, mesh_ima.topology().dim(), 0)
vmap0 = mesh_ima.data().array("parent_vertex_indices", 0)
cmap0 = mesh_ima.data().array("parent_cell_indices", mesh_ima.topology().dim())

        
n = 0
for c in cells(mesh_ima):
  parent_cell = Cell(mesh, cmap0[c.index()])
  surface_marker_ima.array()[c.index()] = cf.array()[parent_cell.index()]
  for f in facets(parent_cell):
    for g in facets(c):
      g_vertices = vmap0[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        boundary_marker_ima.array()[g.index()] = mf.array()[f.index()]
    n=n+1

ds_ima = Measure("ds", domain=mesh_ima, subdomain_data=boundary_marker_ima) 
dx_ima= Measure("dx", domain=mesh_ima, subdomain_data=surface_marker_ima)


#------------------------------------------------------SOLVING HYPERELASTICITY PROBLEM ON THE OUTER 3 LAYERS-----------------------------------------------------------------------------

V = VectorFunctionSpace(mesh_ima, "Lagrange", 1) 
n = FacetNormal(mesh_ima) 

#------------------------------------------------------DEFINING FUNCTIONS NEEDED FOR THE VARIATIONAL FORM--------------------------------------------------------------------------------

du = TrialFunction(V)            # Incremental displacement
u  = Function(V)                 # Deformation vector
v  = TestFunction(V)             # Test function 

#-------------------------------------------------------------SETTING THE KINEMATICS-----------------------------------------------------------------------------------------------------
d = mesh_ima.geometry().dim()
I = Identity(d)             
F = I + grad(u)             
g_int = Identity(d)  
g_med = Identity(d)  
g_ad = Identity(d)  
C_int = (F*inv(g_int)).T*(F*inv(g_int))                   
C_med = (F*inv(g_med)).T*(F*inv(g_med))                   
C_ad = (F*inv(g_ad)).T*(F*inv(g_ad))                   
M = cofac(F)
P = Constant(5.0)
  

Ic_int = tr(C_int)     
Ic_med = tr(C_med)     
Ic_ad = tr(C_ad)     
J_int  = det(F*inv(g_int))
J_med = det(F*inv(g_med))
J_ad = det(F*inv(g_ad))
J_e = det(F)
C_lumen = F.T*F
Ic_lumen = tr(C_lumen)

# Elasticity parameters 
nu = 0.497   
mu1, mu2, mu3 = 27.9, 1.27, 7.56 
eta_int, eta_med, eta_ad = 263.66, 21.6, 38.57
beta_int, beta_med, beta_ad = 170.88, 8.21, 85.03
rho_int, rho_med, rho_ad = 0.51, 0.25, 0.55

# Stored strain energy density 
psi1 = (mu1/2)*(Ic_int - 3) + ((nu)*(mu1)/(1-2*nu))*(J_int-1)**2 - mu1*(ln(J_int)) + (eta_int/beta_int)*(exp(beta_int*((1-rho_int)*(Ic_int - 3)**2))-1)
psi2 = (mu2/2)*(Ic_med - 3) + ((nu)*(mu2)/(1-2*nu))*(J_med-1)**2 - mu2*(ln(J_med)) + (eta_med/beta_med)*(exp(beta_med*((1-rho_med)*(Ic_med - 3)**2))-1)
psi3 = (mu3/2)*(Ic_ad - 3) + ((nu)*(mu3)/(1-2*nu))*(J_ad-1)**2 - mu3*(ln(J_ad)) + (eta_ad/beta_ad)*(exp(beta_ad*((1-rho_ad)*(Ic_ad - 3)**2))-1)

  
Pi = (det(g_int))*(psi1*dx_ima(5)) + (det(g_med))*(psi2*dx_ima(6)) + (det(g_ad))*(psi3*dx_ima(7)) - (0.5)*dot(M*P*(-n), u)*ds_ima(1)

# Compute first variation of Pi (directional derivative about u in the direction of v)
F1 = derivative(Pi, u, v)

# Compute Jacobian of F
Ju = derivative(F1, u, du)

bcs = []

# Solve variational problem
solve(F1 == 0, u, bcs,
  form_compiler_parameters=ffc_options)
  
num_dofs_per_component = int(V.dim()/V.num_sub_spaces())
num_sub_spaces = V.num_sub_spaces()
vector = np.zeros((num_sub_spaces, num_dofs_per_component))
for i in range(num_sub_spaces):
    vector[i] = u.sub(i, deepcopy=True).vector().get_local()
x = V.sub(0).collapse().tabulate_dof_coordinates()
vector = vector.T
print(vector,vector.shape)  # This works perfectly

ALE.move(mesh_ima,u)

#*****************************************Restricting the solution to boundary with the physical tag 1**********************************
 
v_2_d = vertex_to_dof_map(V_orig)
dofs1 = []

for facet in facets(mesh):
    if mf[facet.index()] == 1:  # The boundary to the innermost mesh has physical tag 1
        vertices = facet.entities(0)
        for vertex in vertices:
            dofs1.append(v_2_d[vertex])
            
unique_dofs1 = np.array(list(set(dofs1)), dtype=np.int32)
boundary_coords1 = V_orig.tabulate_dof_coordinates()[unique_dofs1]
u_e = u.vector()[unique_dofs1]
print(u_e, u_e.shape)

#*****************************************Printing the solution at the boundary in a vector form****************************************

num_dofs_per_component_1 = int(V_orig.dim()/V_orig.num_sub_spaces())
num_sub_spaces_1 = V_orig.num_sub_spaces()
vector1 = np.zeros((num_sub_spaces_1, num_dofs_per_component_1))
for i in range(num_sub_spaces_1):
    vector1[i] = u_e.sub(i, deepcopy=True).vector().get_local() # This produces the error
x = V_orig.sub(0).collapse().tabulate_dof_coordinates()
vector1 = vector1.T
print(vector1,vector1.shape)

Is there a problem with the way I restricted the solution to the boundary using a vertex_to_dof_map? Or is the problem somewhere else?

4

The problem is the line above, u_e is a numpy array after assignment, and Therefore does not have the attribute sub.

For any further questions, please make a minimal example of what you want to do, and what goes wrong, as your current example is not following any of the guidelines in Read before posting: How do I get my question answered?