-
Notifications
You must be signed in to change notification settings - Fork 6
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #8 from NGSolve/uzerbinati/fenicsx
Added FEniCSx support
- Loading branch information
Showing
7 changed files
with
209 additions
and
9 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,81 @@ | ||
| ''' | ||
| This module contains all the functions related to wrapping NGSolve meshes to FEniCSx | ||
| We adopt the same docstring conventiona as the FEniCSx project, since this part of | ||
| the package will only be used in combination with FEniCSx. | ||
| ''' | ||
| import typing | ||
| import dolfinx | ||
| import numpy as np | ||
|
|
||
| from mpi4py import MPI as _MPI | ||
|
|
||
| from ngsPETSc import MeshMapping | ||
|
|
||
| # Map from Netgen cell type (integer tuple) to GMSH cell type | ||
| _ngs_to_cells = {(2,3): 2, (2,4):3, (3,4): 4} | ||
|
|
||
| class GeometricModel: | ||
| """ | ||
| This class is used to wrap a Netgen geometric model to a DOLFINx mesh. | ||
| Args: | ||
| geo: The Netgen model | ||
| comm: The MPI communicator to use for mesh creation | ||
| """ | ||
| def __init__(self,geo, comm: _MPI.Comm): | ||
| self.geo = geo | ||
| self.comm = comm | ||
|
|
||
| def model_to_mesh(self, hmax: float, gdim: int = 2, | ||
| partitioner: typing.Callable[ | ||
| [_MPI.Comm, int, int, dolfinx.cpp.graph.AdjacencyList_int32], | ||
| dolfinx.cpp.graph.AdjacencyList_int32] = | ||
| dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.none), | ||
| transform: typing.Any = None, routine: typing.Any = None) -> typing.Tuple[dolfinx.mesh.Mesh, | ||
| dolfinx.cpp.mesh.MeshTags_int32,dolfinx.cpp.mesh.MeshTags_int32]: | ||
| """Given a NetGen model, take all physical entities of the highest | ||
| topological dimension and create the corresponding DOLFINx mesh. | ||
| This function only works in serial, at the moment. | ||
| Args: | ||
| hmax: The maximum diameter of the elements in the triangulation | ||
| model: The NetGen model | ||
| gdim: Geometrical dimension of the mesh | ||
| partitioner: Function that computes the parallel | ||
| distribution of cells across MPI ranks | ||
| transform: PETSc DMPLEX Transformation to be applied to the mesh | ||
| routine: Function to be applied to the mesh after generation | ||
| takes as plan the mesh and the NetGen model and returns the | ||
| same objects after the routine has been applied. | ||
| Returns: | ||
| A DOLFINx mesh for the given NetGen model. | ||
| """ | ||
| # First we generate a mesh | ||
| ngmesh = self.geo.GenerateMesh(maxh=hmax) | ||
| # Apply any ngs routine post meshing | ||
| if routine is not None: | ||
| ngmesh, self.geo = routine(ngmesh, self.geo) | ||
| # Applying any PETSc Transform | ||
| if transform is not None: | ||
| meshMap = MeshMapping(ngmesh) | ||
| transform.setDM(meshMap.plex) | ||
| transform.setUp() | ||
| newplex = transform.apply(meshMap.plex) | ||
| meshMap = MeshMapping(newplex) | ||
| ngmesh = meshMap.ngmesh | ||
| # We extract topology and geometry | ||
| if ngmesh.dim == 2: | ||
| V = ngmesh.Coordinates() | ||
| T = ngmesh.Elements2D().NumPy()["nodes"] | ||
| T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1 | ||
| elif ngmesh.dim == 3: | ||
| V = ngmesh.Coordinates() | ||
| T = ngmesh.Elements3D().NumPy()["nodes"] | ||
| T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1 | ||
| ufl_domain = dolfinx.io.gmshio.ufl_mesh(_ngs_to_cells[(gdim,T.shape[1])],gdim) | ||
| cell_perm = dolfinx.cpp.io.perm_gmsh(dolfinx.cpp.mesh.to_type(str(ufl_domain.ufl_cell())), | ||
| T.shape[1]) | ||
| T = T[:, cell_perm] | ||
| mesh = dolfinx.mesh.create_mesh(self.comm, T, V, ufl_domain, partitioner) | ||
| return mesh |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,68 @@ | ||
| ''' | ||
| This module test the utils.fenicsx class | ||
| ''' | ||
| import pytest | ||
|
|
||
| def test_square_netgen(): | ||
| ''' | ||
| Testing FEniCSx interface with Netgen generating a square mesh | ||
| ''' | ||
| try: | ||
| from mpi4py import MPI | ||
| import ngsPETSc.utils.fenicsx as ngfx | ||
| from dolfinx.io import XDMFFile | ||
| except ImportError: | ||
| pytest.skip(msg="DOLFINx unavailable, skipping FENICSx test") | ||
|
|
||
| from netgen.geom2d import SplineGeometry | ||
| geo = SplineGeometry() | ||
| geo.AddRectangle((0,0),(1,1)) | ||
| geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD) | ||
| domain =geoModel.model_to_mesh(hmax=0.1) | ||
| with XDMFFile(domain.comm, "XDMF/mesh.xdmf", "w") as xdmf: | ||
| xdmf.write_mesh(domain) | ||
|
|
||
| def test_poisson_netgen(): | ||
| ''' | ||
| Testing FEniCSx interface with Netgen generating a square mesh | ||
| ''' | ||
| try: | ||
| import numpy as np | ||
| import ufl | ||
| from dolfinx import fem, mesh | ||
| from dolfinx.fem.petsc import LinearProblem | ||
| from ufl import dx, grad, inner | ||
| from mpi4py import MPI | ||
| from petsc4py.PETSc import ScalarType | ||
| import ngsPETSc.utils.fenicsx as ngfx | ||
| except ImportError: | ||
| pytest.skip(msg="DOLFINx unavailable, skipping FENICSx test") | ||
|
|
||
| from netgen.geom2d import SplineGeometry | ||
| geo = SplineGeometry() | ||
| geo.AddRectangle((0,0),(np.pi,np.pi)) | ||
| geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD) | ||
| msh =geoModel.model_to_mesh(hmax=0.1) | ||
| V = fem.FunctionSpace(msh, ("Lagrange", 2)) | ||
| facetsLR = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1), | ||
| marker=lambda x: np.logical_or(np.isclose(x[0], 0.0), | ||
| np.isclose(x[0], np.pi))) | ||
| facetsTB = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1), | ||
| marker=lambda x: np.logical_or(np.isclose(x[1], 0.0), | ||
| np.isclose(x[1], np.pi))) | ||
| facets = np.append(facetsLR,facetsTB) | ||
| dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets) | ||
| bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V) | ||
| u = ufl.TrialFunction(V) | ||
| v = ufl.TestFunction(V) | ||
| x = ufl.SpatialCoordinate(msh) | ||
| f = ufl.exp(ufl.sin(x[0])*ufl.sin(x[1])) | ||
| a = inner(grad(u), grad(v)) * dx | ||
| L = inner(f, v) * dx | ||
| problem = LinearProblem(a, L, bcs=[bc], | ||
| petsc_options={"ksp_type": "cg", "pc_type": "qr"}) | ||
| problem.solve() | ||
|
|
||
| if __name__ == "__main__": | ||
| test_square_netgen() | ||
| test_poisson_netgen() |