Hi,
I need to use “Real type element” within the mixed function space in dolfinx, as shown below, but I run into the error as attached below, Is there a way I can fix this? or what is the replacement for real type element in dolfinx?. Can anybody help me out with this if you know.
import dolfinx
import numpy as np
from mpi4py import MPI
import dolfinx.fem as fem
from petsc4py import PETSc
from dolfinx.fem.petsc import (apply_lifting, apply_lifting_nest,
assemble_matrix, assemble_matrix_block,
assemble_matrix_nest, assemble_vector,
assemble_vector_block, assemble_vector_nest,
create_matrix, create_matrix_block,
create_matrix_nest, create_vector,
create_vector_block, create_vector_nest, set_bc,
set_bc_nest)
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction,
TestFunction, div, ds, grad, inner,dot, nabla_grad, nabla_div,Identity,dot, dx)
from dolfinx.mesh import CellType, create_box, locate_entities_boundary
from dolfinx.mesh import locate_entities
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, VectorFunctionSpace,
locate_dofs_topological,assemble_scalar,form)
import ufl
from ufl import (Form, SpatialCoordinate, VectorElement, FiniteElement, TensorElement, MixedElement,as_tensor,as_vector,dot,TestFunction, TrialFunction,TestFunctions, TrialFunctions,dx, ds, grad, FacetNormal,inner, max_value,nabla_grad, nabla_div, Identity)
from scipy.optimize import root
from petsc4py.PETSc import ScalarType
import dolfinx.io
from petsc4py import PETSc
from dolfinx.fem import (apply_lifting, locate_dofs_geometrical, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
N = 1
L, B, H = 1,1,1
mesh = dolfinx.mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, B, H])], [N,N,N])
#Create mesh
U = VectorElement("CG", mesh.ufl_cell(), 1)# # displacement vector element
V = FiniteElement("CG", mesh.ufl_cell(), 1) # voltage finite element
V0 = FiniteElement("Real", mesh.ufl_cell(), 0) # Lagrange multiplier
U1, V1 = FunctionSpace(mesh, U), FunctionSpace(mesh, V)
W = FunctionSpace(mesh, MixedElement([U, V, V0]))
num_dofs_global = W.dofmap.index_map.size_global * W.dofmap.index_map_bs
print(f"Number of dofs global: {num_dofs_global}")
# Parameters of Anisotropic Piezomaterial (PIC-255)
c_11 = (1.23e11)
c_12 = (7.67e10)
c_13 = (7.023e10)
c_33 = (9.711e10)
c_44 = (2.226e10)
# Relative permitivity matrix parameters
eps_0 = (8.854e-12) # relative permitivity of free space
eps_11 = (1649 * eps_0)
eps_33 = (1750 * eps_0)
# Coupling matrix parameters
e_13 = (-7.15)
e_33 = (13.75)
e_15 = (11.91)
# Order of matrices similar to ANSYS notation: x, y, z, xy, yz,xz
# IEEE notation: x, y, z, yz, xz, xy
# Elasticity tensor c^E
#IEEE Notation
C = as_tensor([[c_11, c_12, c_13, 0., 0., 0.],
[c_12, c_11, c_13, 0., 0., 0.],
[c_13, c_13, c_33, 0., 0., 0.],
[0., 0., 0., c_44, 0., 0],
[0., 0., 0., 0., c_44, 0.],
[0., 0., 0., 0., 0., (c_11-c_12)/2]])
# piezoelectric coupling tensor e
et_piezo = as_tensor([[0.0, 0.0, 0.0, 0.0, e_15, 0.0],
[0.0, 0.0, 0.0, e_15, 0.0, 0.0],
[e_13, e_13, e_33, 0.0, 0.0, 0.0]])
# transpose form of piezo tensor
e_piezo = as_tensor([[0.0, 0.0, e_13],
[0.0, 0.0, e_13],
[0.0, 0.0, e_33],
[0.0, e_15, 0.0],
[e_15, 0.0, 0.0],
[0.0, 0.0, 0.0]])
# Permittivitats tensor epsilon^S
eps_s = as_tensor([[eps_11, 0., 0.],
[0., eps_11, 0.],
[0., 0., eps_33]])
# rewrite the tensor into a vector using Voigt notation
def strain3voigt(ten):
# FEniCS does not know anything about Voigt notation,
# so one need to access the components directly as eps[0, 0] etc.
return as_vector([ten[0,0],ten[1,1],ten[2,2],2*ten[1,2],2*ten[0,2],2*ten[0,1]])
# rewrite a vector into a tensor using Voigt notation
def voigt3stress(vec):
return as_tensor([[vec[0], vec[5], vec[4]],
[vec[5], vec[1], vec[3]],
[vec[4], vec[3], vec[2]]])
# first define the function for the strain tensor in fenics notation
def B(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#define the electric field vector
def E_field_vector(Ei):
return as_vector([Ei[0], Ei[1], Ei[2]])
#define the electric field in relation to potential
def E_field(v):
return -grad(v)
def Sigma(u):
return voigt3stress(dot(C, strain3voigt(B(u))))
#definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_piezo.E)
def sigma(u, v):
return voigt3stress(dot(C, strain3voigt(B(u))) - e_piezo * E_field_vector(E_field(v)))
# the electric displacements vector
def disp_D(Di):
return as_vector([Di[0], Di[1], Di[2]])
# definition of the different components of my electric displacements equation (D=e^s:Bu+eps^S.E)
# elect_disp_D for coupled linear elasicity and electrostatics
def D_field(u, v):
D = et_piezo * strain3voigt(B(u)) + eps_s * E_field(v)
return disp_D(D)
# Displacement Dirichlet BC
fdim = mesh.topology.dim - 1
def clamped_boundary(x):
return np.isclose(x[1], 0)
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)
u_D = Function(U1)
with u_D.vector.localForm() as loc:
loc.set(0)
bc1 = dirichletbc(u_D, locate_dofs_topological((W.sub(0), U1), fdim, boundary_facets), W.sub(0))
# Voltage Dirichlet BC
def Ground(x):
return np.isclose(x[2], 1)
ground_facets = locate_entities_boundary(mesh, fdim, Ground)
v_D = Function(V1)
with v_D.vector.localForm() as loc:
loc.set(0)
bc2 = dirichletbc(v_D, locate_dofs_topological((W.sub(1),V1), fdim, ground_facets), W.sub(1))
bcs =[bc1,bc2]
# Trial and Test functions
(u, v, v0) = TrialFunctions(W)
(wu, wv, wv0) = TestFunctions(W)
# Boundary Marking
boundaries = [(2, lambda x: np.isclose(x[2], 0))] # coupling surface
# loop through all the boundary conditions to identify the facets
facet_indices, facet_markers = [], []
for (marker, locator) in boundaries:
facets = locate_entities_boundary(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure('ds', domain=mesh, subdomain_data=facet_tag)
n = FacetNormal(mesh)
T = Constant(mesh, ScalarType((0,0,0)))
Tv = Constant(mesh, ScalarType((1e-6))) # surface charge
he = 2*Circumradius(mesh)
alpha = (100)/he
# Weak Formulation
Kuuv = (inner(sigma(u, v), B(wu)))*dx
Kvvu = (inner(-D_field(u, v), grad(wv)))*dx \
+ inner(D_field(u, v), n)* wv * ds(2) + inner(D_field(wu, wv), n)*(v)*ds(2) - alpha*inner((v), wv)*ds(2) - inner(D_field(u, v), n)* wv0 * ds(2)
soln = Function(W)
a = Kuuv + Kvvu
L = inner(T, wu)*ds - inner(Tv, wv)* ds(2) + inner(D_field(wu, wv), n)*(5)*ds(2) + alpha*inner((5), wv)*ds(2)
# Set up the PDE
problem = fem.petsc.LinearProblem(a, L, bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
soln = problem.solve()
Ur = soln.sub(0).collapse()
Vr = soln.sub(1).collapse()
V0 = soln.sub(2).collapse()
error message:
Traceback (most recent call last):
File "/root/3D_PiezoCube/piezocube_FP_Nitsche.py", line 46, in <module>
W = FunctionSpace(mesh, MixedElement([U, V, V0]))
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py", line 448, in __init__
(self._ufcx_element, self._ufcx_dofmap), module, code = jit.ffcx_jit(
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/jit.py", line 206, in ffcx_jit
r = ffcx.codegeneration.jit.compile_elements([ufl_object], parameters=p_ffcx, **p_jit)
File "/usr/local/lib/python3.9/dist-packages/ffcx/codegeneration/jit.py", line 126, in compile_elements
impl = _compile_objects(decl, elements, names, module_name, p, cache_dir,
File "/usr/local/lib/python3.9/dist-packages/ffcx/codegeneration/jit.py", line 232, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, parameters=parameters)
File "/usr/local/lib/python3.9/dist-packages/ffcx/compiler.py", line 103, in compile_ufl_objects
ir = compute_ir(analysis, object_names, prefix, parameters, visualise)
File "/usr/local/lib/python3.9/dist-packages/ffcx/ir/representation.py", line 85, in compute_ir
ir_elements = [
File "/usr/local/lib/python3.9/dist-packages/ffcx/ir/representation.py", line 86, in <listcomp>
_compute_element_ir(e, analysis.element_numbers, finite_element_names)
File "/usr/local/lib/python3.9/dist-packages/ffcx/ir/representation.py", line 122, in _compute_element_ir
basix_element = create_element(ufl_element)
File "/usr/local/lib/python3.9/dist-packages/ffcx/element_interface.py", line 66, in create_element
family_type = basix.finite_element.string_to_family(family_name, element.cell().cellname())
File "/usr/local/lib/python3.9/dist-packages/basix/finite_element.py", line 87, in string_to_family
raise ValueError(f"Unknown element family: {family} with cell type {cell}")
ValueError: Unknown element family: Real with cell type tetrahedron