adding Dirichlet Boundary Normal for disk

doc/fenicsx/doflinx_fem.ipynb

@@ -100,7 +100,7 @@
   },
   "language_info": {
    "name": "python",
-   "version": "3.9.6 (default, Oct 18 2022, 12:41:40) \n[Clang 14.0.0 (clang-1400.0.29.202)]"
+   "version": "3.9.6"
   },
   "vscode": {
    "interpreter": {

src/bvpy/boundary_conditions/boundary.py

@@ -33,7 +33,7 @@ _XYZ = {'x': Symbol('x[0]'),
         'y': Symbol('x[1]'),
         'z': Symbol('x[2]'),
         'all': Symbol('on_boundary'),
-        'near': Function('near')}
+        'near': Function('numpy.isclose')}
 
 
 class Boundary(object):
@@ -84,7 +84,7 @@ class Boundary(object):
 
         # TODO: change this name in to an adequate one
         _x = symbols("x")
-        self._domain = lambdify(_x, self._expr )
+        self._domain = lambdify(_x, self._expr, 'numpy')
 
     def __call__(self):
         """Fetches the fenics object defining the boundary.

src/bvpy/boundary_conditions/dirichlet.py

@@ -97,10 +97,9 @@ class ConstantDirichlet(object):
             #                      method=self._method)
             # boundary_dofs = fem.locate_dofs_geometrical(functionSpace, self._boundary())
             _domain = functionSpace.mesh
-            tdim = _domain.topology.dim
-            _domain.topology.create_connectivity(tdim-1, tdim)
-            boundary_factes = mesh.exterior_facet_indices(_domain.topology)
-            boundary_dofs = fem.locate_dofs_topological(functionSpace, tdim- 1, boundary_factes)
+            fdim = _domain.topology.dim - 1
+            boundary_factes = mesh.locate_entities_boundary(_domain, fdim, self._boundary())
+            boundary_dofs = fem.locate_dofs_topological(functionSpace, fdim, boundary_factes)
             dir = fem.dirichletbc(ScalarType(self._val), boundary_dofs, functionSpace)
             return dir
         else:
@@ -376,7 +375,11 @@ class NormalDirichlet(object):
                + str(self._boundary._string) + ", value: " + str(self._val)+">"
 
     def apply(self, functionSpace):
-        return None
+
+        _fdim = functionSpace.mesh.topology.dim - 1
+        boundary_facets = mesh.locate_entities_boundary(functionSpace.mesh, _fdim, self._boundary())
+        boundary_dofs = fem.locate_dofs_topological(functionSpace, _fdim, boundary_facets)
+        return fem.dirichletbc(boundary_normal(functionSpace, boundary_facets, self._val),boundary_dofs)
         # return fe.DirichletBC(functionSpace,
         #                       boundary_normal(functionSpace.mesh(),
         #                                       scale=self._val),

src/bvpy/domains/geometry.py

@@ -28,9 +28,12 @@ import ufl
 from dolfinx import fem
 from dolfinx.mesh import Mesh
 from dolfinx.mesh import exterior_facet_indices
+import dolfinx as dx
+from petsc4py import PETSc
+from mpi4py import MPI
 
 from bvpy.domains.abstract import AbstractDomain
-
+from bvpy.utils.decorators.boundary_condition import facet_normal
 
 # class Init_orientation_3D(fe.UserExpression):
 class Init_orientation_3D(object):
@@ -172,6 +175,58 @@ def normal(domain, scale=1, interior=None, degree=1, name='normal'):
     # u.name = 'label'
     return u
 
+@facet_normal
+def complex_domain_facet_normal(V, mesh_tag: dx.mesh.meshtags, mesh_tag_id: int, scale: int = 0):
+    """
+    """
+    n = ufl.FacetNormal(V.mesh)
+    u_h = fem.Function(V)
+    u = ufl.TrialFunction(V)
+    v = ufl.TestFunction(V)
+    ds = ufl.ds(domain=V.mesh, subdomain_data=mesh_tag, subdomain_id=mesh_tag_id)
+    a = ufl.inner(u, v) * ds
+    L = ufl.inner(n, v) * ds
+
+    imap = V.dofmap.index_map
+    all_blocks = np.arange(imap.size_local, dtype=np.int32)
+    top_blocks = fem.locate_dofs_topological(V, V.mesh.topology.dim - 1, mesh_tag.find(mesh_tag_id))
+    deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]
+
+    bilinear_form = fem.form(a)
+    pattern = fem.create_sparsity_pattern(bilinear_form)
+    pattern.insert_diagonal(deac_blocks)
+    pattern.assemble()
+    u_0 = fem.Function(V)
+    u_0.vector.set(0.)
+
+    bc_deac = fem.dirichletbc(u_0, deac_blocks)
+    A = dx.cpp.la.petsc.create_matrix(V.mesh.comm, pattern)
+    A.zeroEntries()
+
+    form_coeffs = dx.cpp.fem.pack_coefficients(bilinear_form)
+    form_consts = dx.cpp.fem.pack_constants(bilinear_form)
+    dx.cpp.fem.petsc.assemble_matrix(A, bilinear_form, form_consts, form_coeffs, [bc_deac])
+    if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
+        A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
+        A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
+        dx.cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [bc_deac], 1.0)
+    A.assemble()
+
+    linear_form = fem.form(L)
+    b = fem.petsc.assemble_vector(linear_form)
+    
+    fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
+    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
+    fem.petsc.set_bc(b, [bc_deac])
+
+    solver = PETSc.KSP().create(MPI.COMM_WORLD)
+    solver.setType(PETSc.KSP.Type.BCGS)
+    solver.setOperators(A)
+    solver.solve(b, u_h.vector)
+    u_h.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
+    return u_h
+
+
 
 def get_boundary_facets(mesh):
     """Selects the boundary elements of a given mesh.
@@ -229,7 +284,7 @@ def get_boundary_vertex(mesh):
     return bnd_vertex
 
 
-def boundary_normal(domain, scale=1, name='boundary_normal'):
+def boundary_normal(domain, boundary,  scale=0, name='boundary_normal'):
     """Computes along the boundary the normals within the tangent plane.
 
     Parameters
@@ -287,8 +342,9 @@ def boundary_normal(domain, scale=1, name='boundary_normal'):
 
     # u.rename(name, 'label')
 
-    V = fem.FunctionSpace(mesh, ("CG", 1))
-    boundary_facets = exterior_facet_indices(mesh.topology)
-    u = fem.locate_dofs_topological(V, mesh.topology-1, boundary_facets)
-
+    fdim = domain.mesh.topology.dim - 1
+    mesh_tag_id = 1
+    marked_values = np.full_like(boundary, mesh_tag_id)
+    facet_tag = dx.mesh.meshtags(domain.mesh, fdim, boundary, marked_values)
+    u = complex_domain_facet_normal(domain, facet_tag, mesh_tag_id, scale)
     return u

src/bvpy/utils/decorators/boundary_condition.py

+import functools
+from dolfinx import fem
+
+def facet_normal(fn):
+    @functools.wraps(fn)
+    def func(*args, **kwargs):
+        u_normal = fem.Function(args[0])
+        u_normal.x.array.real[:] = args[3]*fn(*args).x.array.real[:]
+        return u_normal
+    return func
\ No newline at end of file

src/bvpy/utils/visu.py

@@ -239,7 +239,7 @@ Puis generer la figure pyvista.
 Difficultes: Adapter les domaines primitives, pour avoir les proprietes
 topology, cell_type et geometry
 """
-def plot(obj):
+def plot(obj, norm=False):
     """
     plot with pyvista
     """
@@ -262,9 +262,8 @@ def plot(obj):
         logging.warning(f"This type of data {obj} can't be plotted !")
 
     _grid = pyvista.UnstructuredGrid(_topology, _cell_types, _geometry)
-
-    if mesh_u_data:
-        _grid.point_data["u"] = obj.x.array.real
+    if mesh_u_data:        
+        _grid.point_data["u"] = obj.x.array.reshape(_geometry.shape[0], len(obj))
         _grid.set_active_scalars("u")
     elif bdry_u_data:
         _grid.point_data["u"] = obj.g.x.array.real
@@ -272,6 +271,11 @@ def plot(obj):
 
     _plotter = pyvista.Plotter()
     _plotter.add_mesh(_grid, show_edges=True)
+
+    if(norm):
+        glyphs = _grid.glyph(orient="u", factor=20)
+        _plotter.add_mesh(glyphs)        
+
     _plotter.view_xy()
     _plotter.show() 
 

src/bvpy/vforms/elasticity.py

@@ -32,6 +32,8 @@ from ufl import grad, sym, tr
 from bvpy.vforms.abstract import AbstractVform, Element
 from bvpy.utils.tensor import applyElementwise
 from bvpy.domains.geometry import normal
+from petsc4py.PETSc import ScalarType
+
 
 __all__ = ['ElasticForm', 'LinearElasticForm', 'HyperElasticForm']
 
@@ -197,23 +199,26 @@ class LinearElasticForm(ElasticForm):
         The computed stress.
 
         """
-        d = u.geometric_dimension()
+        d = u.ufl_shape[0]
         lmbd = self._parameters['lmbda']
         mu = self._parameters['mu']
         fe_stress = lmbd * tr(self._strain(u)) * ufl.Identity(d) +\
                     2 * mu * self._strain(u)
 
-        if self._plane_stress and d==3:
-            # Use instead FacetNormal from ufl
-            # n = ufl.FacetNormal(mesh)
-            nrl = normal(u.ufl_function_space().mesh)
-            tgt_proj = ufl.Identity(d) - ufl.outer(nrl, nrl)
+        # if self._plane_stress and d==3:
+        #     # Use instead FacetNormal from ufl
+        #     # n = ufl.FacetNormal(mesh)
+        #     nrl = normal(u.ufl_function_space().mesh)
+        #     tgt_proj = ufl.Identity(d) - ufl.outer(nrl, nrl)
 
-            return ufl.dot(tgt_proj, ufl.dot(fe_stress, tgt_proj))
+        #     return ufl.dot(tgt_proj, ufl.dot(fe_stress, tgt_proj))
 
-        else:
-            return fe_stress
+        # else:
+        #     return fe_stress
 
+        # fe_stress = lmbd * ufl.nabla_div(u) * ufl.Identity(u.ufl_shape[0]) + 2*mu*ufl.sym(grad(u)) 
+        return fe_stress
+    
     def construct_form(self, u, v, sol):
 
         mesh = v.ufl_function_space().mesh
@@ -270,7 +275,7 @@ class HyperElasticForm(ElasticForm):
             The computed strain.
 
         """
-        d = u.geometric_dimension()
+        d = u.ufl_shape[0]
         Id = ufl.Identity(d)
         F = Id + ufl.nabla_grad(u)
         C = F.T * F
@@ -308,7 +313,7 @@ class HyperElasticForm(ElasticForm):
             The computed stress.
 
         """
-        d = u.geometric_dimension()
+        d = u.ufl_shape[0]
 
         if model == 'st_venant':
             return self._parameters['lmbda'] * ufl.tr(self._strain(u))\
@@ -328,7 +333,7 @@ class HyperElasticForm(ElasticForm):
             W += self._parameters['mu'] * ufl.tr(E * E)
 
         elif self.model == 'neo_hookean':
-            d = u.geometric_dimension()
+            d = u.ufl_shape[0]
             Id = ufl.Identity(d)
             F = Id + ufl.nabla_grad(sol)
             C = F.T * F