Hello,
I am trying to solve an electrostatic problem in 3D. So first I create the mesh using the opencascade kernel, then I solve for a scalar field (electric potential) and then I apply grad() to the electric potential to obtain a vector field with the electric field.
As the final step, I would like to integrate over a surface to find the amount of charge inside a volume enclosed by the surface.
I receive the following error:
Discontinuous type Jacobian must be restricted.
ERROR:UFL:Discontinuous type Jacobian must be restricted.
The surface I use for the integral is the same I use as a boundary condition for the electric potential. It is basically just a long box.
I checked with the gmsh gui that the surface is closed and all of them have the correct tag. Otherwise I am not sure what could cause this error?
I woul very much appreciate any help or suggestions.
Here is the full code, the part where I try to integrate over a surface is at the bottom.
# https://jsdokken.com/src/tutorial_gmsh.html
# https://fenicsproject.org/olddocs/dolfinx/dev/python/demos/gmsh/demo_gmsh.py.html
mesh_filename = '01_mesh.msh'
resolution = 2 # Normal resolution
resolution_hq = 0.1 # High Quality resolution
# Pick some tags for the elements
SUBS_VOLUME_TAG = 1
AIR_VOLUME_TAG = 2
METAL_WALLS_TAG = 1
MICROSTRIP_WALLS_TAG = 2
SUBS_AIR_WALLS_TAG = 3
# The substrate width, length and thickness
subs_w, subs_h, subs_th = 30, 30, 0.254
# The substrate dielectic constant
subs_er = 2.33
# The air box above width, length and height
air_w, air_h, air_th = subs_w, subs_h, 5
# The microstrip width, length and thickness
metal_w, metal_L, metal_th = 0.7, 20, 0.02
import gmsh
gmsh.initialize()
#%% The first part - prepare the model
# Create a new (and only model)
gmsh.model.add("dut")
# Preapre the substrate and the air volumes
substrate = gmsh.model.occ.addBox(0, -subs_h/2, -subs_th, subs_w, subs_h, subs_th)
air = gmsh.model.occ.addBox(0, -air_h/2, 0, air_w, air_h, air_th)
# Create the metal volume
metal_offset = (subs_w-metal_L)/2
metal = gmsh.model.occ.addBox(metal_offset, -metal_w/2, 0, metal_L, metal_w, metal_th)
# Remove the metal volume from the air volume
# We want to first remove the metal form the air, and then from the substrate, so in the first occ.cut there is removeTool=False
cut_result, _ = gmsh.model.occ.cut([(3, air)], [(3, metal)], removeTool=False)
air = cut_result[0][1]
cut_result, _ = gmsh.model.occ.cut([(3, substrate)], [(3, metal)])
substrate = cut_result[0][1]
# The most important thing, do the boolean fragment - then the mesh is conformal
gmsh.model.occ.fragment([(3, air), (3, substrate)], [ ])
#%% The second part - assign tags
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
# We will try to determine which volume is which and mark them
# later we use this information to assign the dielectric
# The substrate is where z<0, the air where z>0
for v in volumes:
com = gmsh.model.occ.getCenterOfMass(v[0],v[1])
if com[2] < 0:
gmsh.model.addPhysicalGroup(v[0], [v[1]], SUBS_VOLUME_TAG)
gmsh.model.setPhysicalName(v[0], SUBS_VOLUME_TAG, "Substrate")
if com[2]>metal_th:
gmsh.model.addPhysicalGroup(v[0], [v[1]], AIR_VOLUME_TAG)
gmsh.model.setPhysicalName(v[0], AIR_VOLUME_TAG, "Air")
# Now we will try to assign appropariate boundary conditions to the surfaces
surfaces = gmsh.model.occ.getEntities(dim=2)
metal_walls = []
microstrip_walls = []
subs_air_interface = []
for surface in surfaces:
com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
bb = gmsh.model.occ.getBoundingBox(2, surface[1])
# All the surfaces Center-of-Mass with z<0 or z>metal_th are metal box
if com[2] < -subs_th/4 or com[2]>metal_th:
metal_walls.append(surface[1])
elif com[2] > 0 and com[2]<1.25*metal_th:
microstrip_walls.append(surface[1])
# There is one metal surface at z=0, the bottom of the microstirp
elif com[2] > -subs_th/4 and com[2] < metal_th/4 and bb[1]>-metal_w and bb[4]<metal_w:
microstrip_walls.append(surface[1])
# And all the other surfaces at z=0 are
elif com[2] > -subs_th/4 and com[2] < metal_th/4 and bb[1]<-metal_w and bb[4]>metal_w:
subs_air_interface.append(surface[1])
else:
raise Warning("Surfaces with not assigned tag found: ", str(surface))
gmsh.model.addPhysicalGroup(2, metal_walls, METAL_WALLS_TAG)
gmsh.model.setPhysicalName(2, METAL_WALLS_TAG, "PEC Box")
gmsh.model.addPhysicalGroup(2, microstrip_walls, MICROSTRIP_WALLS_TAG)
gmsh.model.setPhysicalName(2, MICROSTRIP_WALLS_TAG, "PEC Microstrip")
gmsh.model.addPhysicalGroup(2, subs_air_interface, SUBS_AIR_WALLS_TAG)
gmsh.model.setPhysicalName(2, SUBS_AIR_WALLS_TAG, "Substrate-Air Interface")
#%% Now let's set the meshing
# First we define from which object the distance matters (from the microstrip)
distance = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance, "FacesList", microstrip_walls)
# So the finest mesh will be up to half of metal_w and then gradually change
# to the regular mesh at the distance of 2x metal_w
threshold = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold, "IField", distance)
gmsh.model.mesh.field.setNumber(threshold, "LcMin", resolution_hq)
gmsh.model.mesh.field.setNumber(threshold, "LcMax", resolution)
gmsh.model.mesh.field.setNumber(threshold, "DistMin", resolution_hq*3)
gmsh.model.mesh.field.setNumber(threshold, "DistMax", resolution_hq*10)
# Set the mesh
gmsh.model.mesh.field.setAsBackgroundMesh(threshold)
#%% Generate the mesh and save it
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)
gmsh.write(mesh_filename)
#%% Load the mesh to DOLFINx
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI
import numpy as np
model_rank = 0
mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim=3)
gmsh.finalize()
# Assign dielectric constant parameters
from dolfinx.fem import (dirichletbc, Expression, Function, FunctionSpace,
VectorFunctionSpace, locate_dofs_topological, Constant)
from dolfinx.fem.petsc import LinearProblem
from ufl import TestFunction, TrialFunction, dot, dx, grad
from petsc4py.PETSc import ScalarType
# The function space, Discrete Galerkin, to assing material properties
Q = FunctionSpace(mesh, ("DG", 0))
# The dielectric constant is a scalar in the whole domain
er = Function(Q)
# Find all substrate cells and assign them er=subs_er
cells = cell_tags.find(SUBS_VOLUME_TAG)
er.x.array[cells] = np.full_like(cells, subs_er, dtype=ScalarType)
# and for the air er=1.0
cells = cell_tags.find(AIR_VOLUME_TAG)
er.x.array[cells] = np.full_like(cells, 1.0, dtype=ScalarType)
# The function space is Continuous Galerkin, because the problem is continuous :D
V = FunctionSpace(mesh, ("CG", 1))
# The dimension of the problem, tdim=3 in this case
tdim = mesh.topology.dim
assert(tdim==3)
# Let's find all the surfaces (or actually facets that compose them) that are assigned to the PEC box, microstrip and the interface between air and dielectric
pec_walls_facets = facet_tags.find(METAL_WALLS_TAG)
pec_microstrip_facets = facet_tags.find(MICROSTRIP_WALLS_TAG)
inteface_facets = facet_tags.find(SUBS_AIR_WALLS_TAG)
# Each of these facets found before has DoFs (degrees of freedom)
# In this case it is points at which a value is assigned, so we want to find all of them so we will be able to fix their value
pec_walls_dofs = locate_dofs_topological(V, 2, pec_walls_facets)
pec_microstrip_dofs = locate_dofs_topological(V, 2, pec_microstrip_facets)
# We create a list of Boundary Conditions (Dirichlet in this case) and assign them electric potential equal to 0 and 1 for the PEC box and microstrip respectievly
bcs = [dirichletbc(ScalarType(0), pec_walls_dofs, V),
dirichletbc(ScalarType(1), pec_microstrip_dofs, V)]
# Preparing the weak formulation
u = TrialFunction(V)
v = TestFunction(V)
# The Poisson equation to find the electric potential
a = er * dot(grad(u), grad(v)) * dx
L = Constant(mesh, ScalarType(0)) * v * dx
# And solving the problem
A = Function(V)
problem = LinearProblem(a, L, u=A, bcs=bcs)
problem.solve()
# Electric field can be found as curl of the electric potential, so first we prepare an appropriate function space
# and then calculate E=grad(A)
W = VectorFunctionSpace(mesh, ("CG", 1))
E = Function(W)
E_expr = Expression(grad(A),
W.element.interpolation_points())
E.interpolate(E_expr)
#%% The final part - calculate the capacitance
# By Gauss's law we have to integrate the electric field over a surface enclosing a volume to find the charge in such volume
# We will take the surface of the microstrip as the surface enclosing the microstrip
# not sure if it is a clever idea, lets find out
from ufl import Measure, FacetNormal
# Create measure
dS = Measure('dS', domain=mesh, subdomain_data=pec_microstrip_facets)
n0 = FacetNormal(mesh)
from dolfinx.fem import assemble_scalar, form
# Integrate
gauss_law = dot(E, n0) * dS
charge = assemble_scalar( form(gauss_law) )
print(charge)
