MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pml.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "pml.hpp"
13
14namespace mfem
15{
16
18 : mesh(mesh_), length(length_)
19{
20 dim = mesh->Dimension();
21 SetBoundaries();
22}
23
24void CartesianPML::SetBoundaries()
25{
26 comp_dom_bdr.SetSize(dim, 2);
27 dom_bdr.SetSize(dim, 2);
28 // initialize
29 for (int i = 0; i < dim; i++)
30 {
31 dom_bdr(i, 0) = infinity();
32 dom_bdr(i, 1) = -infinity();
33 }
34
35 for (int i = 0; i < mesh->GetNBE(); i++)
36 {
37 Array<int> bdr_vertices;
38 mesh->GetBdrElementVertices(i, bdr_vertices);
39 for (int j = 0; j < bdr_vertices.Size(); j++)
40 {
41 for (int k = 0; k < dim; k++)
42 {
43 dom_bdr(k, 0) = std::min(dom_bdr(k, 0), mesh->GetVertex(bdr_vertices[j])[k]);
44 dom_bdr(k, 1) = std::max(dom_bdr(k, 1), mesh->GetVertex(bdr_vertices[j])[k]);
45 }
46 }
47 }
48
49#ifdef MFEM_USE_MPI
50 ParMesh * pmesh = dynamic_cast<ParMesh *>(mesh);
51 if (pmesh)
52 {
53 for (int d=0; d<dim; d++)
54 {
55 MPI_Allreduce(MPI_IN_PLACE, &dom_bdr(d,0), 1,
56 MPITypeMap<real_t>::mpi_type, MPI_MIN,pmesh->GetComm());
57 MPI_Allreduce(MPI_IN_PLACE, &dom_bdr(d,1), 1,
58 MPITypeMap<real_t>::mpi_type, MPI_MAX, pmesh->GetComm());
59 }
60 }
61#endif
62
63 for (int i = 0; i < dim; i++)
64 {
65 comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
66 comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
67 }
68}
69
71 Array<int> * attrPML)
72{
73 int nrelem = mesh_->GetNE();
74 elems.SetSize(nrelem);
75
76 for (int i = 0; i < nrelem; ++i)
77 {
78 elems[i] = 1;
79 bool in_pml = false;
80 Element *el = mesh_->GetElement(i);
81 Array<int> vertices;
82 // Initialize Attribute
83 el->SetAttribute(1);
84 el->GetVertices(vertices);
85 int nrvert = vertices.Size();
86 // Check if any vertex is in the pml
87 for (int iv = 0; iv < nrvert; ++iv)
88 {
89 int vert_idx = vertices[iv];
90 real_t *coords = mesh_->GetVertex(vert_idx);
91 for (int comp = 0; comp < dim; ++comp)
92 {
93 if (coords[comp] > comp_dom_bdr(comp, 1) ||
94 coords[comp] < comp_dom_bdr(comp, 0))
95 {
96 in_pml = true;
97 break;
98 }
99 }
100 }
101 if (in_pml)
102 {
103 elems[i] = 0;
104 el->SetAttribute(2);
105 }
106 }
107 mesh_->SetAttributes();
108
109 if (mesh_->attributes.Size())
110 {
111 if (attrNonPML)
112 {
113 attrNonPML->SetSize(mesh_->attributes.Max());
114 *attrNonPML = 0; (*attrNonPML)[0] = 1;
115
116 }
117 if (attrPML)
118 {
119 attrPML->SetSize(mesh_->attributes.Max());
120 *attrPML = 0;
121 if (mesh_->attributes.Max()>1)
122 {
123 (*attrPML)[1]=1;
124 }
125 }
126 }
127
128}
129
131 std::vector<std::complex<real_t>> &dxs)
132{
133 std::complex<real_t> zi = std::complex<real_t>(0., 1.);
134
135 real_t n = 2.0;
136 real_t c = 5.0;
137 real_t coeff;
138 real_t k = omega * sqrt(epsilon * mu);
139 // Stretch in each direction independently
140 for (int i = 0; i < dim; ++i)
141 {
142 dxs[i] = 1.0;
143 if (x(i) >= comp_dom_bdr(i, 1))
144 {
145 coeff = n * c / k / pow(length(i, 1), n);
146 dxs[i] = real_t(1.0) + zi * real_t(coeff * std::abs(pow(x(i) - comp_dom_bdr(i,
147 1), n - 1.0)));
148 }
149 if (x(i) <= comp_dom_bdr(i, 0))
150 {
151 coeff = n * c / k / pow(length(i, 0), n);
152 dxs[i] = real_t(1.0) + zi * real_t(coeff * std::abs(pow(x(i) - comp_dom_bdr(i,
153 0), n - 1.0)));
154 }
155 }
156}
157
158// acoustics UW PML coefficients functions
159// |J|
161{
162 int dim = pml->dim;
163 std::vector<std::complex<real_t>> dxs(dim);
164 std::complex<real_t> det(1.0,0.0);
165 pml->StretchFunction(x, dxs);
166 for (int i=0; i<dim; ++i) { det *= dxs[i]; }
167 return det.real();
168}
169
171{
172 int dim = pml->dim;
173 std::vector<std::complex<real_t>> dxs(dim);
174 std::complex<real_t> det(1.0,0.0);
175 pml->StretchFunction(x, dxs);
176 for (int i=0; i<dim; ++i) { det *= dxs[i]; }
177 return det.imag();
178}
179
181{
182 int dim = pml->dim;
183 std::vector<std::complex<real_t>> dxs(dim);
184 std::complex<real_t> det(1.0,0.0);
185 pml->StretchFunction(x, dxs);
186 for (int i=0; i<dim; ++i) { det *= dxs[i]; }
187 return det.imag()*det.imag() + det.real()*det.real();
188}
189
190// J^T J / |J|
192 DenseMatrix & M)
193{
194 int dim = pml->dim;
195 std::vector<std::complex<real_t>> dxs(dim);
196 std::complex<real_t> det(1.0,0.0);
197 pml->StretchFunction(x, dxs);
198 for (int i = 0; i<dim; ++i) { det *= dxs[i]; }
199
200 M=0.0;
201 for (int i = 0; i<dim; ++i)
202 {
203 M(i,i) = (pow(dxs[i], real_t(2))/det).real();
204 }
205}
206
208 DenseMatrix & M)
209{
210 int dim = pml->dim;
211 std::vector<std::complex<real_t>> dxs(dim);
212 std::complex<real_t> det = 1.0;
213 pml->StretchFunction(x, dxs);
214 for (int i = 0; i<dim; ++i) { det *= dxs[i]; }
215
216 M=0.0;
217 for (int i = 0; i<dim; ++i)
218 {
219 M(i,i) = (pow(dxs[i], real_t(2))/det).imag();
220 }
221}
222
224 DenseMatrix & M)
225{
226 int dim = pml->dim;
227 std::vector<std::complex<real_t>> dxs(dim);
228 std::complex<real_t> det = 1.0;
229 pml->StretchFunction(x, dxs);
230 for (int i = 0; i<dim; ++i) { det *= dxs[i]; }
231
232 M=0.0;
233 for (int i = 0; i<dim; ++i)
234 {
235 std::complex<real_t> a = pow(dxs[i], real_t(2))/det;
236 M(i,i) = a.imag() * a.imag() + a.real() * a.real();
237 }
238}
239
240
241// Maxwell PML coefficients
243 DenseMatrix &M)
244{
245 int dim = pml->dim;
246 std::vector<std::complex<real_t>> dxs(dim);
247 std::complex<real_t> det(1.0, 0.0);
248 pml->StretchFunction(x, dxs);
249
250 for (int i = 0; i < dim; ++i) { det *= dxs[i]; }
251
252 M = 0.0;
253 for (int i = 0; i < dim; ++i)
254 {
255 M(i, i) = (det / pow(dxs[i], real_t(2))).real();
256 }
257}
258
260 DenseMatrix &M)
261{
262 int dim = pml->dim;
263 std::vector<std::complex<real_t>> dxs(dim);
264 std::complex<real_t> det = 1.0;
265 pml->StretchFunction(x, dxs);
266
267 for (int i = 0; i < dim; ++i) { det *= dxs[i]; }
268
269 M = 0.0;
270 for (int i = 0; i < dim; ++i)
271 {
272 M(i, i) = (det / pow(dxs[i], real_t(2))).imag();
273 }
274}
275
277 DenseMatrix &M)
278{
279 int dim = pml->dim;
280 std::vector<std::complex<real_t>> dxs(dim);
281 std::complex<real_t> det = 1.0;
282 pml->StretchFunction(x, dxs);
283
284 for (int i = 0; i < dim; ++i) { det *= dxs[i]; }
285
286 M = 0.0;
287 for (int i = 0; i < dim; ++i)
288 {
289 std::complex<real_t> a = det / pow(dxs[i], real_t(2));
290 M(i, i) = a.real()*a.real() + a.imag()*a.imag();
291 }
292}
293
294} // namespace mfem
Dynamic 2D array using row-major layout.
Definition array.hpp:372
void SetSize(int m, int n)
Definition array.hpp:385
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Class for setting up a simple Cartesian PML region.
Definition pml.hpp:19
void SetAttributes(Mesh *mesh_, Array< int > *attrNonPML=nullptr, Array< int > *attrPML=nullptr)
Mark element in the PML region.
Definition pml.cpp:70
void StretchFunction(const Vector &x, std::vector< std::complex< real_t > > &dxs)
PML complex stretching function.
Definition pml.cpp:130
CartesianPML(Mesh *mesh_, const Array2D< real_t > &length_)
Definition pml.cpp:17
real_t epsilon
Definition pml.hpp:49
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Definition element.hpp:58
Mesh data type.
Definition mesh.hpp:56
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition mesh.hpp:1442
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition mesh.cpp:1604
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
real_t a
Definition lissajous.cpp:41
real_t detJ_r_function(const Vector &x, CartesianPML *pml)
PML stretching functions: See https://doi.org/10.1006/jcph.1994.1159.
Definition pml.cpp:160
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:45
void Jt_J_detJinv_r_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:191
void abs_Jt_J_detJinv_2_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:223
void detJ_Jt_J_inv_r_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:242
void Jt_J_detJinv_i_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:207
float real_t
Definition config.hpp:43
real_t abs_detJ_2_function(const Vector &x, CartesianPML *pml)
Definition pml.cpp:180
real_t detJ_i_function(const Vector &x, CartesianPML *pml)
Definition pml.cpp:170
void abs_detJ_Jt_J_inv_2_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:276
void detJ_Jt_J_inv_i_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:259