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.