ok this is the code for the mesh, just launch it to create the file.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 22 15:28:10 2024
"""
import sys
try:
import gmsh
except ImportError:
print("This demo requires gmsh to be installed")
sys.exit(0)
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
# WARNING: unit length is 1.0 nm
c=50## coté carré
esp=20## espacement carré avec 2dmat
ld=600## long domaine
hd=300## haut domaine
mesh_res_r_min = 1
mesh_res_r_max = 20
gmsh.initialize()
gmsh.model.add("syst")
# Geometry
#### le carré à 20nm du 2dmat
p1 = gmsh.model.occ.addPoint(-c/2, esp, 0.0, mesh_res_r_min)
p2 = gmsh.model.occ.addPoint(c/2, esp, 0.0, mesh_res_r_min)
p3 = gmsh.model.occ.addPoint(c/2, esp+c, 0.0, mesh_res_r_min)
p4 = gmsh.model.occ.addPoint(-c/2, esp+c, 0.0, mesh_res_r_min)
#### le domaine entier 600 x 300 nm
p5 = gmsh.model.occ.addPoint(-ld/2, 0.0, 0.0, mesh_res_r_min)
p6 = gmsh.model.occ.addPoint(+ld/2, 0.0, 0.0, mesh_res_r_min)
p7 = gmsh.model.occ.addPoint(+ld/2, hd, 0.0, mesh_res_r_max)
p8 = gmsh.model.occ.addPoint(-ld/2, hd, 0.0, mesh_res_r_max)
### carré
line1 = gmsh.model.occ.addLine(p1, p2)
line2 = gmsh.model.occ.addLine(p2, p3)
line3 = gmsh.model.occ.addLine(p3, p4)
line4 = gmsh.model.occ.addLine(p4, p1)
### domaine
line5 = gmsh.model.occ.addLine(p5, p6)
line6 = gmsh.model.occ.addLine(p6, p7)
line7 = gmsh.model.occ.addLine(p7, p8)
line8 = gmsh.model.occ.addLine(p8, p5)
carre = gmsh.model.occ.addCurveLoop([line1, line2, line3, line4])
surf_carre = gmsh.model.occ.addPlaneSurface([carre])
domaine = gmsh.model.occ.addCurveLoop([line5, line6, line7, line8])
surf_dom = gmsh.model.occ.addPlaneSurface([domaine])
surf_poisson = gmsh.model.occ.fragment([(2, surf_dom)], [(2, surf_carre)], removeTool=False)
# #####
gmsh.model.occ.synchronize()
bottom_line = [gmsh.model.getBoundary(surf_poisson[1][0])[0]] ## après multiple modif
gmsh.model.mesh.setSize(bottom_line, mesh_res_r_min) ### bonne syntax mais ne change pas la resol du mesh
#### /!\ /!\ #### ne marchait pas car il fallait ajouter que les points au coin du bas avait eux aussi cette resolution
##############################################################################
boundary_dirichlet = gmsh.model.addPhysicalGroup(1, [line1, line2,line3,line4])
boundary_neumann_air_haut = gmsh.model.addPhysicalGroup(1, [line7])
boundary_neumann_air_gauche = gmsh.model.addPhysicalGroup(1, [line8])
boundary_neumann_air_droite = gmsh.model.addPhysicalGroup(1, [line6])
boundary_neumann_2dmat = gmsh.model.addPhysicalGroup(1, [line5])
gmsh.model.setPhysicalName(1, boundary_neumann_2dmat, "boundary_2dmat")
##############################################################################
### ajouter Permittivité relative au matériaux
dielectric_poisson = gmsh.model.addPhysicalGroup(2,surf_poisson[0][0]) ##### le [0][0] correspond à la surface 2-1 dom-carré, [1][0] c'est les lignes je crois
gmsh.model.setPhysicalName(2, dielectric_poisson, "dielectric_poisson")
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate()
gmsh.write("mesh2d_liu_test.msh")
and this is the program, i can’t really make it smaller, most important part is from line 200 to 249 : " Solver Poisson".
line 243 is what i want to make work. g = script (value of u at BC)
in this mesh the Neumann BC is at y=0
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu May 23 15:18:00 2024
"""
#### environnement fenew
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import dolfinx.fem.petsc
import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import meshio
from petsc4py import PETSc
###### ajouté depuis https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import vtk_mesh
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
import pyvista
from dolfinx import mesh, fem, io, nls, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import CellType, create_rectangle, locate_entities, locate_entities_boundary
##############################################################################
def fct_dens_chrg(x,u):
u=-13.295*pow(u,2)*np.sign(u)
return u
##############################################################################
class PoissonSolverBase:
"""
Class that defines basic quantities useful for all
Poisson solvers
"""
def __init__(self):
self.proc_id = MPI.COMM_WORLD.rank
self.markers_1d = {}
self.markers_2d = {}
self.markers_3d = {}
self.bcs = []
def __create_meshio(self, mesh, cell_type, prune_z=False):
"""
Select parts of an existing meshio mesh and create an existing meshio mesh
"""
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
def reads_gmsh_mesh(self, name_of_file):
"""
Reads gmsh mesh
"""
if self.proc_id == 0:
# Read in mesh
self.mesh_meshio = meshio.read(name_of_file)
def convert_to_xdmf(self, name_out_file, element_type, prune_z):
"""
Convert gmsh mesh to xdmf
"""
if self.proc_id == 0:
# Create and save one file for the mesh, and one file for the facets
mesh_tmp =self. __create_meshio(self.mesh_meshio, element_type, prune_z=prune_z)
meshio.write(name_out_file, mesh_tmp)
def create_mesh(self, name_of_file):
"""
Read the H5-XDMF files and generate the mesh
"""
with XDMFFile(MPI.COMM_WORLD, name_of_file, "r") as xdmf:
self.mesh_geom = xdmf.read_mesh(name="Grid")
self.mesh_geom_tags = xdmf.read_meshtags(self.mesh_geom, name="Grid")
self.mesh_geom.topology.create_connectivity(self.mesh_geom.topology.dim, self.mesh_geom.topology.dim-1)
def create_mesh_boundaries(self, name_of_file):
"""
Read the meshtags for boundaries
"""
with XDMFFile(MPI.COMM_WORLD, name_of_file, "r") as xdmf:
self.mesh_boundaries_tags = xdmf.read_meshtags(self.mesh_geom, name="Grid")
def set_dielectric_constant(self, eps_val_markers):
"""
Set the dielectric constant function
"""
self.dielectric_cons = fem.Function(self.Q)
for values in eps_val_markers:
cells_tmp = self.mesh_geom_tags.find(values[0])
self.dielectric_cons.x.array[cells_tmp] = np.full_like(cells_tmp, values[1], dtype=ScalarType)
def set_dirichlet_bc(self, marker, value):
"""
Set Dirichlet boundary conditions from marker
"""
facets = self.mesh_boundaries_tags.find(marker)
dofs = fem.locate_dofs_topological(self.V, self.mesh_geom.topology.dim-1, facets)
self.bcs.append(fem.dirichletbc(ScalarType(value), dofs, self.V)) ## je crois c'est ça à changer pour neumann value ce sra la value de g, et on veut g=0
######################## Neumann BC ###########################
def set_neumann_bc(self, markers, value, test_fct):
"""
Set Neumann boundary conditions from marker
"""
if type(markers)==int :
markers=[markers]
facet_indices, facet_markers = [], []
for marker in markers:
facets = self.mesh_boundaries_tags.find(marker)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(self.mesh_geom, self.mesh_geom.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
dss = ufl.Measure("ds", domain=self.mesh_geom, subdomain_id = facet_markers[0], subdomain_data = facet_tag) # avec subdomain_data = facet_tag on a assertion error
Formul = -1*inner(value, test_fct) * dss
return Formul
##############################################################################
def define_measure_on_boundary(self, markers):
facet_indices, facet_markers = [], []
for marker in markers:
facets = self.mesh_boundaries_tags.find(marker)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(self.mesh_geom, self.mesh_geom.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
return ufl.Measure("dS", domain=self.mesh_geom, subdomain_data=facet_tag)
def solve(self):
"""
Solve the linear problem
"""
self.uh = self.problem.solve()
def save_solution_xdmf(self, name_of_file):
"""
Save the solution to xdfm format
"""
with io.XDMFFile(self.mesh_geom.comm, name_of_file, "w") as file:
file.write_mesh(self.mesh_geom)
file.write_function(self.uh)
def visualize_solution(self):
"""
Visualize the solution with pyvista
"""
try:
import pyvista
cells, types, x = plot.vtk_mesh(self.V) ### plot.create_vtk obsolete^
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = self.uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show()
except ModuleNotFoundError:
print("'pyvista' is required to visualise the solution")
print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
###############################################################################
###############################################################################
def fct_for_uh(self):
tdim = self.mesh_geom.topology.dim-2# got the Dofs
cells1 = locate_entities_boundary(self.mesh_geom, tdim, lambda m: m[1] == 0) ## is exactly the indexes where x[1]>=59.9
coords= self.V.tabulate_dof_coordinates()
xx=coords[cells1]# coordinate of these dofs
uu=self.uh.eval(xx,cells1)# value of u at these dofs
uu=fct_dens_chrg(xx[0],uu)## pour pybinding il faudra r soit x[0]
self.uh.vector[cells1]=uu
return self.uh
###############################################################################
################### Solver Poisson #################
class system_to_solve(PoissonSolverBase):
"""
Class that solves the Poisson equation for Liu non linear system
"""
def __init__(self, name_of_gmsh_file, eps, bias):
super().__init__()
self.reads_gmsh_mesh(name_of_gmsh_file)
self.convert_to_xdmf("mesh_system.xdmf", "triangle", True)
self.convert_to_xdmf("mesh_system_boundaries.xdmf", "line", True)
self.create_mesh("mesh_system.xdmf")
self.create_mesh_boundaries("mesh_system_boundaries.xdmf")
# Define markers
self.markers_1d["square"] = 1 # carre
self.markers_1d["left"] = 4 # bord gauche
self.markers_1d["right"] = 2 # bord droite
self.markers_1d["top"] = 3 # bord haut
self.markers_1d["mat2d"] = 5 ## bord bas ici
self.markers_2d["air"] = 6 # le milieu air
# Create function space
self.Q = fem.functionspace(self.mesh_geom, ("DG", 0)) ## c'est pour .set_dielectric_constant
self.V = fem.functionspace(self.mesh_geom, ("CG", 1)) ##
#################### u v ##### pour NEWTON nonlinear :
self.uh = fem.Function(self.V)
self.v = ufl.TestFunction(self.V)
# Create a function dielectric_cons for the dielectric constant with different value in every region (one region here)
eps_markers = [(self.markers_2d["air"], eps[0])]
self.set_dielectric_constant(eps_markers)
########################## F ##############################
F= self.dielectric_cons * inner(grad(self.uh), grad(self.v)) * dx
###Apply Dirichlet boundary conditions where needed using markers
self.set_dirichlet_bc(self.markers_1d["square"], bias)
####### Neumann border ext + 2dmat
# g2d = -13.295*pow(self.uh,2)*ufl.sign(self.uh)## explicit neumann BC ( WORKING)
g2d = self.fct_for_uh() # that's my objective script dependant Neumann BC not working (not far?)
F+=self.set_neumann_bc(self.markers_1d["mat2d"],g2d,self.v)
self.problem = NonlinearProblem(F,self.uh, bcs=self.bcs)
##########################################################################
bias = 1.5 # potential at the square
eps = [1.0] # dielec constant, air
to_solve = system_to_solve("mesh2d_liu_test.msh", eps, bias)
###### non linear Newton
solver = NewtonSolver(MPI.COMM_WORLD, to_solve.problem)
solver.rtol = 1e-6
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}setInitialGuessNonzero"] = "True"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(to_solve.uh)
assert (converged)
print(f"Number of iterations: {n:d}")
########### SAVE AND FIGURE
to_solve.save_solution_xdmf("sol_poisson_liu_system.xdmf")
to_solve.visualize_solution()
cells, types, x = plot.vtk_mesh(to_solve.V)
UH = to_solve.uh.x.array.real
y=x[:,1]
## les x en y = 0
ind=np.where(y<1E-10)[0] # 601 values, i took 1nm step for 2dmat for x = -300 : +300nm, so looks good
## UH for y=0
UHX = UH[ind]
### we take x at y=0
X = x[ind,0]
sort=np.argsort(X) # index in increasing order for X
X=X[sort]### applied to X
UHX=UHX[sort] ### applied to UHX
plt.plot(X,UHX,'-')
plt.title('Potentiel(x) y=0')
hope it helps