I am quit new in FeniCs. I want to know how can I generate polygonal meshes and then define triangle subdomain on each mesh?
I would be thankful if sb answer my question.
1 Like
I would suggest using GMSH, which DOLFINx has an interface to, see for instance
https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_gmsh.html
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html
http://jsdokken.com/FEniCS23-tutorial/src/mesh_generation.html
https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html
Thank you Dokken for your respond.
But in GMSH how can I generate polygon. It only can generate triangle and rectangles?
You can crate polygons in GMSH by creating line loops with line segments, then in turn triangulate these into triangles.
FEniCS does not support polyhedral cells, but triangles, tetrahedrons, quadrilaterals and hexahedrons.
Thank you for your quick answer dokken.
Hi all,
I made a collection of subdomains within my domain. These subdomains have a set of facets that they coincide in two different neighbourhood subdomains but they are not identical. I want to compute the jump the function over these facets. But as the index of the facets are different I don’t know how to do that. That is I want to compute:
u^(-)ds^(-)-u^(+)ds^(+)
Do you have any idea that how can I do it?
What do you mean by this?
Without a minimal reproducible code it is quite hard to fathom what you are trying to do.
When I generated the mesh using gmsh, I employed distinct indices for the points and facets of two adjacent subdomains. First, I generated the polygonal meshes , then refined it into triangles (as shown in the image). The rationale behind this mesh construction lies in my intention to implement a Recovered Discontinuous Galerkin (DG) method on the polygonal mesh, while utilizing a continuous Galerkin (CG) method for the refined mesh (which is recovered by CG).
Since I am using the recovered finite element so necessarily I don’t need to work with a polygonal mesh in FeniCS. But for the DG method I need to compute the jump over facets of the polygonal mesh. My idea was just use the DG method. In Fact by adopting this mesh indices for the subdomains, the assemble of the matrix on the polygonal mesh results in a blocked matrix similar to that in the DG method. But still I need to compute the jump over polygonal facets (subdomain facets) to determine the consistency and penalty terms.
You would define a MeshFunction over your polygonal domains with for instance a single tag indicating the boundary between any polygon. Then, you should send this meshfunction into the dS measure (Interior facet mesh, as subdomain_data) and call the integral restriction.
Imagine you have your MeshFunction ff indicating the polygonal domains, then you would call jump(u)*dS(subdomain_data=ff, subdomain_id=id_of_polygonal_facets).
Note that if the order of “+” and minus matters, you would need to do some extra work, i.e.
Hi Dokken,
Thank you so much for your answer. I wrote my mesh with a unique tag.
Then added the subdomain boundaries to dS. Now, I can compute jump(u)dS. But for the variational form I need to compute
inner(jump(u, n), avg(grad(v))) * dS
and
jump(u) jump(v)) * dS
which does not work. I mean when I multiply it by the jump or average of test function I got the following error
TypeError Traceback (most recent call last)
Cell In[4], line 251
249 a = lhs(F)
250 L = rhs(F)
→ 251 problem = LinearProblem(a, L, bcs=[bc], petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
253 uh = problem.solve()
254 ### Interpolate Exact solution ##################################################################
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py:588, in LinearProblem.init(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
556 def init(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBC] = ,
557 u: typing.Optional[_Function] = None,
558 petsc_options: typing.Optional[dict] = None,
559 form_compiler_options: typing.Optional[dict] = None,
560 jit_options: typing.Optional[dict] = None):
561 “”“Initialize solver for a linear variational problem.
562
563 Args:
(…)
586 “pc_factor_mat_solver_type”: “mumps”})
587 “””
→ 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)
File /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:183, in form.._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 /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:175, in form.._form(form)
171 # Subdomain markers (possibly empty list for some integral types)
172 subdomains = {_ufl_to_dolfinx_domain[key]: get_integration_domains(
173 _ufl_to_dolfinx_domain[key], subdomain_data[0]) for (key, subdomain_data) in sd.get(domain).items()}
→ 175 f = ftype(module.ffi.cast(“uintptr_t”, module.ffi.addressof(ufcx_form)), V, coeffs,
176 constants, subdomains, mesh)
177 return Form(f, ufcx_form, code)
TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. dolfinx.cpp.fem.Form_float64(spaces: List[dolfinx::fem::FunctionSpace], integrals: Dict[dolfinx::fem::IntegralType, List[Tuple[int, object, numpy.ndarray[numpy.int32]]]], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, mesh: dolfinx.cpp.mesh.Mesh_float64 = None)
2. dolfinx.cpp.fem.Form_float64(form: int, spaces: List[dolfinx::fem::FunctionSpace], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], subdomains: Dict[dolfinx::fem::IntegralType, List[Tuple[int, numpy.ndarray[numpy.int32]]]], mesh: dolfinx.cpp.mesh.Mesh_float64)
Invoked with: <cdata ‘uintptr_t’ 140566361245312>, [<dolfinx.cpp.fem.FunctionSpace_float64 object at 0x7fd85473faf0>, <dolfinx.cpp.fem.FunctionSpace_float64 object at 0x7fd85473faf0>], [<dolfinx.cpp.fem.Function_float64 object at 0x7fd8081a93f0>], , {<IntegralType.cell: 0>: , <IntegralType.interior_facet: 2>: 1}, <dolfinx.cpp.mesh.Mesh_float64 object at 0x7fd80905dd50>
Please provide a minimal reproducible example with 3x` format,
Ie
```
# add code here
```
Hi, Dokken, This is a part of my code:
alpha = fem.Constant(msh, default_real_type(6.0 * k**2))
h = CellDiameter(msh)
msh, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()
grid = pyvista.UnstructuredGrid(*vtk_mesh(msh, msh.topology.dim))
num_local_cells = msh.topology.index_map(msh.topology.dim).size_local
grid.cell_data[“Marker”] = ct.values[ct.indices < num_local_cells]
marker_values=grid.set_active_scalars(“Marker”)
dS= Measure(“dS”, domain=msh, subdomain_data=ct)
a = epsilondot(grad(u), grad(v))dx+dot(B,grad(u))vdx - inner(avg(grad(u)), jump(v, n)) * dSinner(jump(u, n), avg(grad(v))) * dS+ (alpha / avg(h)) epsilon* inner(jump(u), jump(v)) * dS*dot(B,n)dot((u(‘+’)-u(‘-’)),outer(v)) dS
FF = f * v * dx
LL = dot(f,v)*dx
F=a-LL
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=[bc], petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
uh = problem.solve(
when I want to compute consistency term and penalty term it makes problem. In fact I only want the average and jump be calculated over the polygonal facets (which are the boundary of the subdomains).
Just giving an unformatted subset of the code is not making it reproducible.
If you arent willing to make an MWE, it is unlikely that people will have time to work through your code.
Hi, Dokken, I am sorry. I misunderstood you. Here you can find my mesh
And this is an example of what I am doing
import sys
import pyvista
from dolfinx.plot import vtk_mesh
from dolfinx.io import XDMFFile, gmshio
from dolfinx.io.gmshio import model_to_mesh
sys.path.insert(0, 'mesh')
import numpy as np
from dolfinx import fem
from dolfinx import default_real_type, mesh
import ufl
from ufl import (FacetNormal, SpatialCoordinate, as_ufl, CellDiameter, FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
div, dot, dS, ds, dx, inner, inner, outer, lhs, grad, nabla_grad, rhs, sym, jump, avg)
from dolfinx.fem import (Constant, Function, FunctionSpace, form, petsc,
assemble_scalar, dirichletbc, form, locate_dofs_topological, locate_dofs_geometrical)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, LinearProblem
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
currentmesh = "mesh"
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.open(currentmesh + ".msh")
model = gmsh.model
model.setCurrent(currentmesh)
print('Model ' + model.getCurrent() + ' (' + str(model.getDimension()) + 'D)')
# Define mesh
gdim = 2
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
model_rank = 0
n = FacetNormal(msh)
k=10
degree=1
V = FunctionSpace(msh, ("CG", degree))
uD = Function(V)
uD.interpolate(lambda x: 0* x[0])
tdim = msh.topology.dim # Highest dimension of the elements (faces)
fdim = tdim - 1 # Degree of the boundary facets (lines)
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(msh.topology) # Find the boundary lines
boundary_dofs = locate_dofs_topological(V=V, entity_dim=fdim, entities=boundary_facets) # Find the boundary degrees of freedom
bc = dirichletbc(uD, boundary_dofs)
u = TrialFunction(V)
v = TestFunction(V)
V = FunctionSpace(msh, ("CG", 1))
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), degree)
V13 = FunctionSpace(msh, P2)
def Velocity(x):
values1 = np.zeros((2, x.shape[1]))
values1[0] = 2+0*x[0]
values1[1] = 1+0*x[0]
return values1
B= Function(V13)
B.interpolate(Velocity)
epsilon=10.0**(-5)
alpha = fem.Constant(msh, default_real_type(6.0 * k**2))
h = CellDiameter(msh)
# Here I add the boundaries of the subdomains to dS
msh, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()
grid = pyvista.UnstructuredGrid(*vtk_mesh(msh, msh.topology.dim))
num_local_cells = msh.topology.index_map(msh.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices < num_local_cells]
marker_values=grid.set_active_scalars("Marker")
dS= Measure("dS", domain=msh, subdomain_data=ct)
a = epsilon*dot(grad(u), grad(v))*dx+dot(B,grad(u))*v*dx - inner(avg(grad(u)), jump(v, n)) * dS- inner(jump(u, n), avg(grad(v))) * dS+ (alpha / avg(h))* epsilon* inner(jump(u), jump(v)) * dS- dot(B,n)*dot((u('+')-u('-')),outer(v))* ds
f= Function(V)
f.interpolate(lambda x: (epsilon*(2*np.pi*np.pi*np.sin(np.pi*x[0])*np.sin(np.pi*x[1]) - 2*np.pi*np.pi*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])*np.cos(x[0]) + np.pi*np.cos(np.pi*x[1])*np.sin(np.pi*x[0])*np.sin(x[0]))+2*np.pi*np.cos(np.pi*x[0])*np.sin(np.pi*x[1])+np.pi*np.sin(np.pi*x[0])*np.cos(np.pi*x[1])+np.sin(np.pi*x[0])*np.sin(np.pi*x[1])))
FF = f * v * dx
LL = dot(f,v)*dx
F=a-LL
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
As I mentioned above I need to compute the jump of some terms such as
inner(avg(grad(u)), jump(v, n)) * dS- inner(jump(u, n), avg(grad(v))) * dS+ (alpha / avg(h))* epsilon* inner(jump(u), jump(v)) * dS
Over the boundary of the subdomains. But I got the above error in the variational form.
ct is not the right subdomain data to provide. ct is a tag for each cell, while the dS measure need a marker of facets (the third argument coming out from mode_to_mesh).
I would also suggest you try to make your example more minimal. For instance, the following code gets past the form compilation
from mpi4py import MPI
import dolfinx.fem.petsc
import ufl
import numpy as np
msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
degree = 1
V = dolfinx.fem.FunctionSpace(msh, ("CG", degree))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
dS= ufl.Measure("dS", domain=msh)
epsilon = dolfinx.fem.Constant(msh, 1.0)
B = ufl.as_vector((1.0, 2.0))
alpha = dolfinx.fem.Constant(msh, 1.0)
n = ufl.FacetNormal(msh)
h = 2*ufl.Circumradius(msh)
a = epsilon*ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx+ufl.dot(B,ufl.grad(u))*v*ufl.dx - ufl.inner(ufl.avg(ufl.grad(u)), ufl.jump(v, n)) * dS- ufl.inner(ufl.jump(u, n), ufl.avg(ufl.grad(v))) * dS+ (alpha / ufl.avg(h))* epsilon* ufl.inner(ufl.jump(u), ufl.jump(v)) * dS- ufl.dot(B,n)*ufl.dot((u('+')-u('-')),ufl.outer(v))* ufl.ds
f= dolfinx.fem.Function(V)
f.interpolate(lambda x: (epsilon*(2*np.pi*np.pi*np.sin(np.pi*x[0])*np.sin(np.pi*x[1]) - 2*np.pi*np.pi*np.cos(np.pi*x[0])*np.cos(np.pi*x[1])*np.cos(x[0]) + np.pi*np.cos(np.pi*x[1])*np.sin(np.pi*x[0])*np.sin(x[0]))+2*np.pi*np.cos(np.pi*x[0])*np.sin(np.pi*x[1])+np.pi*np.sin(np.pi*x[0])*np.cos(np.pi*x[1])+np.sin(np.pi*x[0])*np.sin(np.pi*x[1])))
FF = f * v * ufl.dx
LL = ufl.dot(f,v)*ufl.dx
F=a-LL
a = ufl.lhs(F)
L = ufl.rhs(F)
dolfinx.fem.form(a)
dolfinx.fem.form(L)
Always try to reduce your example to something that can be replicated on a built in mesh, and please post the full error trace formatted with
```bash
Add error trace here
```
Hi Dokken,
This is my error:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[5], line 40
38 a = lhs(F)
39 L = rhs(F)
---> 40 problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
42 uh = problem.solve()
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py:588, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
556 def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBC] = [],
557 u: typing.Optional[_Function] = None,
558 petsc_options: typing.Optional[dict] = None,
559 form_compiler_options: typing.Optional[dict] = None,
560 jit_options: typing.Optional[dict] = None):
561 """Initialize solver for a linear variational problem.
562
563 Args:
(...)
586 "pc_factor_mat_solver_type": "mumps"})
587 """
--> 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)
File /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:133, in form.<locals>._form(form)
131 # Extract subdomain data from UFL form
132 sd = form.subdomain_data()
--> 133 domain, = list(sd.keys()) # Assuming single domain
134 # Check that subdomain data for each integral type is the same
135 for data in sd.get(domain).values():
ValueError: too many values to unpack (expected 1)
You are loading in two different meshes
You should only do this once (“with gdim=gdim”), and use facet_markers to access the boundaries you want to use.
Thank you so much for your answer. It finally worked for me. I have another question I also need to compute
(\boldsymbol{\beta}\cdot \boldsymbol{n}(u_{h}^{+}-u_{h}^{-}),v_{h}^{+})_{\partial_{-} K}
where \beta is the vector field and K is any cell (I want to apply inflow boundary on all cells) and
\partial_{-}K:= \{\ x \in \partial K: \beta(x) \cdot \bf{n}(x) < 0\}.
I used this
flux = ufl.dot(B, n) # The flux you want to integrate
cond = ufl.lt(flux, 0.0) # "lower than" condition, returns True where flux < 0 and False elsewhere
conditional_flux = ufl.conditional(cond, flux, 0.0)
But it only works on the boundary and when I wanted to use it for all the internal facets I got the variational error.
that is when I use it as
conditional_flux*(u('+')-u('-'))*v('+')* dS
I got the error:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[4], line 135
133 a = lhs(F)
134 L = rhs(F)
--> 135 problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
137 uh = problem.solve()
138 #u_PPM = functions1.PPM(V, msh, bcs, boundary_dofs, uh, FF, a, A, B, u, v, n, dS, h_avg, penalty_parameter, jump_grad_v, jump_v, a_max=k, omega=0.03, max_iterations=500, N=N, epsilon=10.0**(-5))[0]
139 #u_PPM = functions1.positivity(u_PPM, k, V)[0]
140 #uh=u_PPM
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py:588, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
556 def __init__(self, a: ufl.Form, L: ufl.Form, bcs: typing.List[DirichletBC] = [],
557 u: typing.Optional[_Function] = None,
558 petsc_options: typing.Optional[dict] = None,
559 form_compiler_options: typing.Optional[dict] = None,
560 jit_options: typing.Optional[dict] = None):
561 """Initialize solver for a linear variational problem.
562
563 Args:
(...)
586 "pc_factor_mat_solver_type": "mumps"})
587 """
--> 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)
File /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-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 /usr/lib/petsc/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 /usr/lib/python3/dist-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 Jacobian must be restricted.
This question is unrelated to the initial problem, and a new posts should be opened wit a minimal reproducible example
Hi Dokken,
Hoverer I added the boundary of the subdomains to the internal facets, but it seems that still it is not computing the jump between two subdomains, Actually as a simple example I used
import sys
import pyvista
import math
from dolfinx.io import XDMFFile, gmshio
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.plot import vtk_mesh
sys.path.insert(0, 'meshes')
#import gmsh
import matplotlib.pyplot as plt
import numpy as np
import numpy.linalg as npla
from scipy.sparse.linalg import spsolve
import scipy.sparse as sps
import ufl
from dolfinx import fem
from dolfinx import default_real_type, mesh
from ufl import FacetNormal, SpatialCoordinate, ds, dx, grad, CellDiameter, inner, dot, as_ufl #Unified form language pieces
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement, TensorElement,
as_vector, div, dot, dS,ds, dx, inner, outer, lhs, grad, nabla_grad, rhs, sym, jump, avg)
from dolfinx.fem import (Expression, Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
import numpy
import gmsh
from mpi4py import MPI
from petsc4py import PETSc
#import functions1
import dolfinx
from dolfinx.io import XDMFFile
#from dolfinx.plot import create_vtk_mesh
import dolfinx.io
from dolfinx import io
########################################
currentmesh = "mesh9" # Change N and the filename for different refinements
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.open(currentmesh + ".msh")
model = gmsh.model
model.setCurrent(currentmesh)
print('Model ' + model.getCurrent() + ' (' + str(model.getDimension()) + 'D)')
# Define mesh
msh, cell_markers, facet_markers = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
k = 1.
n = FacetNormal(msh)
################################################
degree=1
V = FunctionSpace(msh, ("CG", degree))
################################################################
tdim = msh.topology.dim
fdim = tdim - 1
u, v = TrialFunction(V), TestFunction(V)
msh.topology.create_connectivity(fdim, tdim)
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
facet_indices1, facet_markers1 = [], []
fdim = msh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(msh, fdim, locator)
facet_indices1.append(facets)
facet_markers1.append(np.full_like(facets, marker))
facet_indices1 = np.hstack(facet_indices1).astype(np.int32)
facet_markers1 = np.hstack(facet_markers1).astype(np.int32)
sorted_facets = np.argsort(facet_indices1)
facet_tag = meshtags(msh, fdim, facet_indices1[sorted_facets], facet_markers1[sorted_facets])
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
ds = Measure("ds", domain=msh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
u_1 = Function(V)
u_1.interpolate(lambda x: 0*x[1]**0)
boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_1),
BoundaryCondition("Dirichlet", 2, u_1),
BoundaryCondition("Dirichlet", 3, u_1),
BoundaryCondition("Dirichlet", 4, u_1),
]
#####################################################################################################################
jump_grad_u = dot(grad(u)('+'),n('+'))-dot(grad(u)('-'),n('-'))
jump_grad_v = dot(grad(v)('+'),n('+'))-dot(grad(v)('-'),n('-'))
jump_u = u('+')*n('+')+u('-')*n('-')
jump_v = v('+')*n('+')+v('-')*n('-')
epsilon=10.0**(-2)
alpha = dolfinx.fem.Constant(msh, 1.0)
h = CellDiameter(msh)
dS= Measure("dS", domain=msh, subdomain_data=facet_markers)
###################################################################################
a = epsilon*dot(grad(u), grad(v))*dx+ u * v * dx- epsilon*inner(avg(grad(v)), jump_u) * dS + epsilon*inner(jump_v, avg(grad(u))) * dS+epsilon*(alpha / avg(h))* epsilon* inner(jump_u, jump_v) * dS#+epsilon*(alpha / avg(h))* epsilon* inner(jump_u, jump_v) * dS
LL= dolfinx.fem.Constant(msh, 1.0) * v * dx
F=a-LL
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
Here I added the facets of the subdomains as
dS= Measure("dS", domain=msh, subdomain_data=facet_markers)
And this is the answer I got
which does not compute any jump between two subdomain.
I used this mesh