Please see below the clean code below for your reference.
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx import default_scalar_type
from ufl import (FacetNormal, FiniteElement, as_matrix, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint
This is mesh
cl__1 = 0.001;
Point(1) = {0.01, 0.07000000000000001, 0, cl__1};
Point(2) = {0.03, 0.07000000000000001, 0, cl__1};
Point(3) = {0.03, 0.08, 0, cl__1};
Point(4) = {0.01, 0.08, 0, cl__1};
Point(5) = {0.04, 0.07000000000000001, 0, cl__1};
Point(6) = {0.06, 0.07000000000000001, 0, cl__1};
Point(7) = {0.06, 0.08, 0, cl__1};
Point(8) = {0.04, 0.08, 0, cl__1};
Point(9) = {0.07000000000000001, 0.07000000000000001, 0, cl__1};
Point(10) = {0.09, 0.07000000000000001, 0, cl__1};
Point(11) = {0.09, 0.08, 0, cl__1};
Point(12) = {0.07000000000000001, 0.08, 0, cl__1};
Point(13) = {0, 0.06, 0, cl__1};
Point(14) = {0.1, 0.06, 0, cl__1};
Point(15) = {0.1, 0.09, 0, cl__1};
Point(16) = {0, 0.09, 0, cl__1};
Point(17) = {0, 0, 0, cl__1};
Point(18) = {0.1, 0, 0, cl__1};
Point(19) = {0.15, 0.1, 0, cl__1};
Point(20) = {-0.05, 0.1, 0, cl__1};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Line(5) = {8, 5};
Line(6) = {5, 6};
Line(7) = {6, 7};
Line(8) = {7, 8};
Line(9) = {12, 9};
Line(10) = {9, 10};
Line(11) = {10, 11};
Line(12) = {11, 12};
Line(13) = {13, 14};
Line(14) = {14, 15};
Line(15) = {15, 16};
Line(16) = {16, 13};
Line(17) = {17, 18};
Line(18) = {18, 19};
Line(19) = {19, 20};
Line(20) = {20, 17};
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Curve Loop(3) = {9, 10, 11, 12};
Plane Surface(3) = {3};
Curve Loop(4) = {1, 2, 3, 4, -8, -7, -6, -5, -12, -11, -10, -9, -15, -14, -13, -16};
Plane Surface(4) = {4};
Curve Loop(5) = {16, 13, 14, 15, -20, -19, -18, -17};
Plane Surface(5) = {5};
Physical Line("bottom") = {17};
Physical Line("right") = {18};
Physical Line("top") = {19};
Physical Line("left") = {20};
Physical Line("EF") = {13};
Physical Line("FG") = {14};
Physical Line("GH") = {19};
Physical Line("HE") = {16};
Physical Line("IJ") = {2};
Physical Line("JK") = {3};
Physical Line("KL") = {4};
Physical Line("LI") = {1};
Physical Line("MN") = {6};
Physical Line("NO") = {7};
Physical Line("OP") = {8};
Physical Line("PM") = {5};
Physical Line("QR") = {10};
Physical Line("RS") = {11};
Physical Line("ST") = {12};
Physical Line("TQ") = {9};
Physical Surface("my_surface") = {5};
Physical Surface("EFGH") = {4};
Physical Surface("IJKL") = {1};
Physical Surface("MNOP") = {2};
Physical Surface("QRST") = {3};
Field[1] = Box;
Field[1].Thickness = 0;
Field[1].VIn = 0;
Field[1].VOut = 0;
Field[1].XMax = 0;
Field[1].XMin = 0;
Field[1].YMax = 0;
Field[1].YMin = 0;
Field[1].ZMax = 0;
Field[1].ZMin = 0;
Here is the code:
# Mesh
mesh_path = "Mesh_fenics.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(mesh_path, MPI.COMM_WORLD, gdim=2)
# Check for Measures
import pyvista
from dolfinx import plot
pyvista.set_jupyter_backend("static")
topology, cell_types, geometry = plot.vtk_mesh(domain, 2)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
p = pyvista.Plotter()
p.add_mesh(grid,show_edges=True)
p.view_xy()
p.show_axes()
p.show()
from dolfinx.plot import vtk_mesh
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "mtgaurav.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(cell_markers, domain.geometry)
import ufl
gdim = 2
def strain(u, repr ="vectorial"):
eps_t = sym(grad(u))
if repr =="vectorial":
return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
elif repr =="tensorial":
return eps_t
E = fem.Constant(domain, 200000000000.0)
nu = fem.Constant(domain,0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu, 0.0],[0.0, 0.0, mu]])
def stress(u, repr ="vectorial"):
sigv = dot(C, strain(u))
if repr =="vectorial":
return sigv
elif repr =="tensorial":
return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])
# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))
#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx
# Applying a uniform pressure
T = fem.Constant(domain, 1000000.0)
n = FacetNormal(domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_markers)
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_markers)
L_form = dot(T*n,u_) * dS(7)
#Boundary Conditions (simply supported on top and bottom edges)
V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))
ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(uy0, top_dofs, V.sub(1)), fem.dirichletbc(ux0, right_dofs, V.sub(0))]
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0))]
WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
WPsolve.solve()
V0 = fem.functionspace(domain, ("DG", 0, (3,)))
sig_exp = fem.Expression(stress(u), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)
import pyvista
from dolfinx import plot
pyvista.set_jupyter_backend("static")
topology, cell_types, geometry = plot.vtk_mesh(domain, 2)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
p = pyvista.Plotter()
p.add_mesh(grid,show_edges=True)
p.view_xy()
p.show_axes()
p.show()
vtk = io.VTKFile(domain.comm, "wpex5.pvd", "u")
vtk.write_function(u, 0)
vtk.write_function(sig, 0)
vtk.close()
Here is the error outcome
--------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[42], line 68
64 #bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(uy0, top_dofs, V.sub(1)), fem.dirichletbc(ux0, right_dofs, V.sub(0))]
65 #bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
66 bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0))]
---> 68 WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
69 WPsolve.solve()
72 V0 = fem.functionspace(domain, ("DG", 0, (3,)))
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/petsc.py:590, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
588 self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
589 self._A = create_matrix(self._a)
--> 590 self._L = _create_form(L, form_compiler_options=form_compiler_options, jit_options=jit_options)
591 self._b = create_vector(self._L)
593 if u is None:
594 # Extract function space from TrialFunction (which is at the
595 # end of the argument list as it is numbered as 1, while the
596 # Test function is numbered as 0)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/forms.py:188, in form(form, dtype, form_compiler_options, jit_options)
185 return list(map(lambda sub_form: _create_form(sub_form), form))
186 return form
--> 188 return _create_form(form)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
181 return form argument"""
182 if isinstance(form, ufl.Form):
--> 183 return _form(form)
184 elif isinstance(form, collections.abc.Iterable):
185 return list(map(lambda sub_form: _create_form(sub_form), form))
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/fem/forms.py:141, in form.<locals>._form(form)
139 if mesh is None:
140 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 141 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
142 form_compiler_options=form_compiler_options,
143 jit_options=jit_options)
145 # For each argument in form extract its function space
146 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
51 @functools.wraps(local_jit)
52 def mpi_jit(comm, *args, **kwargs):
53
54 # Just call JIT compiler when running in serial
55 if comm.size == 1:
---> 56 return local_jit(*args, **kwargs)
58 # Default status (0 == ok, 1 == fail)
59 status = 0
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
202 # Switch on type and compile, returning cffi object
203 if isinstance(ufl_object, ufl.Form):
--> 204 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
205 elif isinstance(ufl_object, ufl.FiniteElementBase):
206 r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:199, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
197 except Exception:
198 pass
--> 199 raise e
201 obj, module = _load_objects(cache_dir, module_name, form_names)
202 return obj, module, (decl, impl)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:190, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
187 for name in form_names:
188 decl += form_template.format(name=name)
--> 190 impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
191 cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
192 except Exception as e:
193 try:
194 # remove c file so that it will not timeout next time
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/codegeneration/jit.py:260, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
256 import ffcx.compiler
258 # JIT uses module_name as prefix, which is needed to make names of all struct/function
259 # unique across modules
--> 260 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
262 ffibuilder = cffi.FFI()
263 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
264 extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/compiler.py:97, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
95 # Stage 1: analysis
96 cpu_time = time()
---> 97 analysis = analyze_ufl_objects(ufl_objects, options)
98 _print_timing(1, time() - cpu_time)
100 # Stage 2: intermediate representation
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/analysis.py:88, in analyze_ufl_objects(ufl_objects, options)
85 else:
86 raise TypeError("UFL objects not recognised.")
---> 88 form_data = tuple(_analyze_form(form, options) for form in forms)
89 for data in form_data:
90 elements += [convert_element(e) for e in data.unique_sub_elements]
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/analysis.py:88, in <genexpr>(.0)
85 else:
86 raise TypeError("UFL objects not recognised.")
---> 88 form_data = tuple(_analyze_form(form, options) for form in forms)
89 for data in form_data:
90 elements += [convert_element(e) for e in data.unique_sub_elements]
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ffcx/analysis.py:163, in _analyze_form(form, options)
160 complex_mode = "_Complex" in options["scalar_type"]
162 # Compute form metadata
--> 163 form_data = ufl.algorithms.compute_form_data(
164 form,
165 do_apply_function_pullbacks=True,
166 do_apply_integral_scaling=True,
167 do_apply_geometry_lowering=True,
168 preserve_geometry_types=(ufl.classes.Jacobian,),
169 do_apply_restrictions=True,
170 do_append_everywhere_integrals=False, # do not add dx integrals to dx(i) in UFL
171 complex_mode=complex_mode)
173 # If form contains a quadrature element, use the custom quadrature scheme
174 custom_q = None
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py:328, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
326 # Propagate restrictions to terminals
327 if do_apply_restrictions:
--> 328 form = apply_restrictions(form)
330 # If in real mode, remove any complex nodes introduced during form processing.
331 if not complex_mode:
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/apply_restrictions.py:180, in apply_restrictions(expression)
177 integral_types = [k for k in integral_type_to_measure_name.keys()
178 if k.startswith("interior_facet")]
179 rules = RestrictionPropagator()
--> 180 return map_integrand_dags(rules, expression,
181 only_integral_type=integral_types)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:73, in map_integrand_dags(function, form, only_integral_type, compress)
71 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
72 """Map integrand dags."""
---> 73 return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
74 form, only_integral_type)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:25, in map_integrands(function, form, only_integral_type)
23 """Apply transform(expression) to each integrand expression in form, or to form if it is an Expr."""
24 if isinstance(form, Form):
---> 25 mapped_integrals = [map_integrands(function, itg, only_integral_type)
26 for itg in form.integrals()]
27 nonzero_integrals = [itg for itg in mapped_integrals
28 if not isinstance(itg.integrand(), Zero)]
29 return Form(nonzero_integrals)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:25, in <listcomp>(.0)
23 """Apply transform(expression) to each integrand expression in form, or to form if it is an Expr."""
24 if isinstance(form, Form):
---> 25 mapped_integrals = [map_integrands(function, itg, only_integral_type)
26 for itg in form.integrals()]
27 nonzero_integrals = [itg for itg in mapped_integrals
28 if not isinstance(itg.integrand(), Zero)]
29 return Form(nonzero_integrals)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:33, in map_integrands(function, form, only_integral_type)
31 itg = form
32 if (only_integral_type is None) or (itg.integral_type() in only_integral_type):
---> 33 return itg.reconstruct(function(itg.integrand()))
34 else:
35 return itg
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py:73, in map_integrand_dags.<locals>.<lambda>(expr)
71 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
72 """Map integrand dags."""
---> 73 return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
74 form, only_integral_type)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/corealg/map_dag.py:34, in map_expr_dag(function, expression, compress, vcache, rcache)
15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
16 """Apply a function to each subexpression node in an expression DAG.
17
18 If the same funtion is called multiple times in a transformation
(...)
32 The result of the final function call
33 """
---> 34 result, = map_expr_dags(function, [expression], compress=compress,
35 vcache=vcache,
36 rcache=rcache)
37 return result
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/corealg/map_dag.py:98, in map_expr_dags(function, expressions, compress, vcache, rcache)
96 # Cache miss: Get transformed operands, then apply transformation
97 if cutoff_types[v._ufl_typecode_]:
---> 98 r = handlers[v._ufl_typecode_](v)
99 else:
100 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/apply_restrictions.py:104, in RestrictionPropagator.reference_value(self, o)
102 f, = o.ufl_operands
103 assert f._ufl_is_terminal_
--> 104 g = self(f)
105 if isinstance(g, Restricted):
106 side = g.side()
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/corealg/multifunction.py:95, in MultiFunction.__call__(self, o, *args)
93 def __call__(self, o, *args):
94 """Delegate to handler function based on typecode of first argument."""
---> 95 return self._handlers[o._ufl_typecode_](o, *args)
File ~/anaconda3/envs/bluemira-dolmpc/lib/python3.11/site-packages/ufl/algorithms/apply_restrictions.py:59, in RestrictionPropagator._require_restriction(self, o)
57 """Restrict a discontinuous quantity to current side, require a side to be set."""
58 if self.current_restriction is None:
---> 59 raise ValueError(f"Discontinuous type {o._ufl_class_.__name__} must be restricted.")
60 return o(self.current_restriction)
ValueError: Discontinuous type Argument must be restricted.