From 414a9f3a28d5f998e9dd3bcbfc8a46ce8cb3ca3a Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 28 May 2024 11:58:21 +0100 Subject: [PATCH] Revert "Added PML support" This reverts commit 9107826e5b115a225aecabbe7e0281216d7d5b1d. --- ngsPETSc/__init__.py | 3 +- ngsPETSc/utils/firedrake/hierarchies.py | 3 +- ngsPETSc/utils/firedrake/pml.py | 78 ------------------------- 3 files changed, 2 insertions(+), 82 deletions(-) delete mode 100644 ngsPETSc/utils/firedrake/pml.py diff --git a/ngsPETSc/__init__.py b/ngsPETSc/__init__.py index 97edab7..8ecdcd5 100644 --- a/ngsPETSc/__init__.py +++ b/ngsPETSc/__init__.py @@ -16,8 +16,7 @@ if firedrake: from ngsPETSc.utils.firedrake.meshes import * from ngsPETSc.utils.firedrake.hierarchies import * - from ngsPETSc.utils.firedrake.pml import * - __all__ = __all__ + ["FiredrakeMesh", "NetgenHierarchy", "PML"] + __all__ = __all__ + ["FiredrakeMesh", "NetgenHierarchy"] #FEniCSx try: diff --git a/ngsPETSc/utils/firedrake/hierarchies.py b/ngsPETSc/utils/firedrake/hierarchies.py index e7c9e28..0d56060 100644 --- a/ngsPETSc/utils/firedrake/hierarchies.py +++ b/ngsPETSc/utils/firedrake/hierarchies.py @@ -1,6 +1,5 @@ ''' -This module contains all the functions related to creating Firedrake mesh hierarchies -from Netgen/NGSolve meshes. +This module contains all the functions related ''' try: import firedrake as fd diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py deleted file mode 100644 index 57ca986..0000000 --- a/ngsPETSc/utils/firedrake/pml.py +++ /dev/null @@ -1,78 +0,0 @@ -''' -This file contains all the functions related to creating a -Perfectly Matched Layer (PML) in Firedrake. -''' - -from firedrake import Constant, Function, solve, DirichletBC, TrialFunction, TestFunction -from firedrake import inner, grad, dx, conditional, real -from firedrake.__future__ import interpolate -from collections.abc import Iterable -class PML: - """ - This class creates a Perfectly Matched Layer (PML) in Firedrake. - :arg mesh: The Firedrake mesh. - :arg V: The function space of the Helmholtz problem. - :arg pmls: A list of touple of the from: - (pml_region_name, pml_inner_boundary_name, pml_outer_boundary_name). - :arg k: The wavenumber of the Helmholtz problem, it can be a Constant, a Function, - or a list, if it is a list we assume each element of the list the wavenumber - for one of the PML regions. - :arg alpha: The attenuation coefficient of the PML, it can be a Constant, a Function, - or a list, if it is a list we assume each element of the list the attenuation - coefficient for one of the PML regions. - :arg solver_parameters: The solver parameters for the PML problem. - """ - def __init__(self, mesh, V, pml_regions, k, alpha=Constant(1j), solver_parameters=None): - """ - This function initializes the PML object. - """ - self.mesh = mesh - self.V = V - self.pml_regions = pml_regions - self.k = k - self.alpha = alpha - if solver_parameters is None: - solver_parameters = {"ksp_type":"preonly", "pc_type":"lu", - "pc_factor_mat_solver_type":"mumps"} - if not hasattr(self.mesh, "netgen_mesh"): - raise ValueError("The mesh must be a Netgen mesh.") - labels1 = dict(map(lambda i,j : (i,j+1) , mesh.netgen_mesh.GetRegionNames(dim=1), - range(len(mesh.netgen_mesh.GetRegionNames(dim=1))))) - labels2 = dict(map(lambda i,j : (i,j+1) , mesh.netgen_mesh.GetRegionNames(dim=2), - range(len(mesh.netgen_mesh.GetRegionNames(dim=2))))) - if not isinstance(k, (Constant, Function, list)): - raise ValueError("The wavenumber must be a Constant, a Function, or a list.") - else: - if not isinstance(k, Iterable): - k = [k]*len(pml_regions) - if not isinstance(alpha, (Constant, Function, list)): - raise ValueError("The attenuation coefficient must be a Constant, a Function, or a list.") - else: - if not isinstance(alpha, Iterable): - alpha = [alpha]*len(pml_regions) - - #Construct the PML for each PML region - self.pmls = [] - for region in self.pml_regions: - if region[0] not in labels2: - raise ValueError("The PML region name must be a valid region name.") - if region[1] not in labels1: - raise ValueError("The PML inner boundary name must be a valid region name.") - if region[2] not in labels1: - raise ValueError("The PML outer boundary name must be a valid region name.") - #Construct the weight function for the PML - u = TrialFunction(V) - v = TestFunction(V) - F = inner(grad(u), grad(v))*dx(labels2[region[0]]) - for i in range(1, len(self.mesh.netgen_mesh.GetRegionNames(dim=2))+1): - if i != labels2[region[0]]: - F += inner(grad(u), grad(v))*dx(i) - L = inner(Constant(1),v)*dx(2) - bcs = [DirichletBC(V, 1, labels1[region[1]]), DirichletBC(V, 0, labels1[region[2]])] - dalet = Function(V) - solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) - sigma = interpolate(1+(alpha/k)*dalet, V) - self.pmls = self.pmls + [conditional(real(sigma) < 1e-12, 1, sigma)] - - -