MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
mesh-fitting.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
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 "mesh-optimizer.hpp"
13#include "../common/mfem-common.hpp"
14
15using namespace std;
16using namespace mfem;
17using namespace common;
18
19// Used for exact surface alignment
21{
22 const int dim = x.Size();
23 if (dim == 2)
24 {
25 const real_t xc = x(0) - 0.5, yc = x(1) - 0.5;
26 const real_t r = sqrt(xc*xc + yc*yc);
27 return r-0.25;
28 }
29 else
30 {
31 const real_t xc = x(0) - 0.5, yc = x(1) - 0.5, zc = x(2) - 0.5;
32 const real_t r = sqrt(xc*xc + yc*yc + zc*zc);
33 return r-0.3;
34 }
35}
36
37real_t in_circle(const Vector &x, const Vector &x_center, real_t radius)
38{
39 Vector x_current = x;
40 x_current -= x_center;
41 real_t dist = x_current.Norml2();
43 {
44 return 1.0;
45 }
46 else if (dist == radius)
47 {
48 return 0.0;
49 }
50
51 return -1.0;
52}
53
55{
56 real_t phi_t = x(1) + (a-b)*x(0)/l - a;
57 return (phi_t <= 0.0) ? 1.0 : -1.0;
58}
59
61{
62 real_t phi_p1 = (x(0)-h-t/2) - k*x(1)*x(1);
63 real_t phi_p2 = (x(0)-h+t/2) - k*x(1)*x(1);
64 return (phi_p1 <= 0.0 && phi_p2 >= 0.0) ? 1.0 : -1.0;
65}
66
68{
69 real_t dx = std::abs(x(0) - xc);
70 real_t dy = std::abs(x(1) - yc);
71 return (dx <= w/2 && dy <= h/2) ? 1.0 : -1.0;
72}
73
74// Fischer-Tropsch like geometry
76{
77 // Circle
78 Vector x_circle1(2);
79 x_circle1(0) = 0.0;
80 x_circle1(1) = 0.0;
81 real_t in_circle1_val = in_circle(x, x_circle1, 0.2);
82
83 real_t r1 = 0.2;
84 real_t r2 = 1.0;
85 real_t in_trapezium_val = in_trapezium(x, 0.05, 0.1, r2-r1);
86
87 real_t return_val = max(in_circle1_val, in_trapezium_val);
88
89 real_t h = 0.4;
90 real_t k = 2;
91 real_t t = 0.15;
92 real_t in_parabola_val = in_parabola(x, h, k, t);
93 return_val = max(return_val, in_parabola_val);
94
95 real_t in_rectangle_val = in_rectangle(x, 0.99, 0.0, 0.12, 0.35);
96 return_val = max(return_val, in_rectangle_val);
97
98 real_t in_rectangle_val2 = in_rectangle(x, 0.99, 0.5, 0.12, 0.28);
99 return_val = max(return_val, in_rectangle_val2);
100 return return_val;
101}
102
103real_t in_cube(const Vector &x, real_t xc, real_t yc, real_t zc, real_t lx,
104 real_t ly, real_t lz)
105{
106 real_t dx = std::abs(x(0) - xc);
107 real_t dy = std::abs(x(1) - yc);
108 real_t dz = std::abs(x(2) - zc);
109 return (dx <= lx/2 && dy <= ly/2 && dz <= lz/2) ? 1.0 : -1.0;
110}
111
112real_t in_pipe(const Vector &x, int pipedir, Vector x_pipe_center,
113 real_t radius, real_t minv, real_t maxv)
114{
115 Vector x_pipe_copy = x_pipe_center;
116 x_pipe_copy -= x;
117 x_pipe_copy(pipedir-1) = 0.0;
118 real_t dist = x_pipe_copy.Norml2();
119 real_t xv = x(pipedir-1);
120 if (dist < radius && xv > minv && xv < maxv)
121 {
122 return 1.0;
123 }
124 else if (dist == radius || (xv == minv && dist < radius) || (xv == maxv &&
126 {
127 return 0.0;
128 }
129
130 return -1.0;
131}
132
134{
135 return r1 + r2 - std::pow(r1*r1 + r2*r2, 0.5);
136}
137
139{
140 return r1 + r2 + std::pow(r1*r1 + r2*r2, 0.5);
141}
142
144{
145 return r_intersect(r1, -r2);
146}
147
149{
150 Vector xcc(3);
151 xcc = 0.5;
152 real_t cube_x = 0.25*2;
153 real_t cube_y = 0.25*2;
154 real_t cube_z = 0.25*2;
155
156 real_t in_cube_val = in_cube(x, xcc(0), xcc(1), xcc(2), cube_x, cube_y, cube_z);
157
158 Vector x_circle_c(3);
159 x_circle_c = 0.5;
160
162 real_t in_sphere_val = in_circle(x, x_circle_c, sphere_radius);
163 real_t in_return_val = std::min(in_cube_val, in_sphere_val);
164
165 int pipedir = 1;
166 Vector x_pipe_center(3);
167 x_pipe_center = 0.5;
171 real_t in_pipe_x = in_pipe(x, pipedir, x_pipe_center, pipe_radius, xmin, xmax);
172
173 in_return_val = std::min(in_return_val, -1*in_pipe_x);
174
175 pipedir = 2;
176 in_pipe_x = in_pipe(x, pipedir, x_pipe_center, pipe_radius, xmin, xmax);
177 in_return_val = std::min(in_return_val, -1*in_pipe_x);
178
179 pipedir = 3;
180 in_pipe_x = in_pipe(x, pipedir, x_pipe_center, pipe_radius, xmin, xmax);
181 in_return_val = std::min(in_return_val, -1*in_pipe_x);
182
183 return in_return_val;
184}
185
186#ifdef MFEM_USE_MPI
188{
189 const int dim = pmesh->Dimension();
190 for (int i = 0; i < pmesh->GetNBE(); i++)
191 {
192 mfem::Array<int> dofs;
193 pmesh->GetNodalFESpace()->GetBdrElementDofs(i, dofs);
194 mfem::Vector bdr_xy_data;
195 mfem::Vector dof_xyz(dim);
196 mfem::Vector dof_xyz_compare;
197 mfem::Array<int> xyz_check(dim);
198 for (int j = 0; j < dofs.Size(); j++)
199 {
200 for (int d = 0; d < dim; d++)
201 {
202 dof_xyz(d) = x(pmesh->GetNodalFESpace()->DofToVDof(dofs[j], d));
203 }
204 if (j == 0)
205 {
206 dof_xyz_compare = dof_xyz;
207 xyz_check = 1;
208 }
209 else
210 {
211 for (int d = 0; d < dim; d++)
212 {
213 if (std::abs(dof_xyz(d)-dof_xyz_compare(d)) < 1.e-10)
214 {
215 xyz_check[d] += 1;
216 }
217 }
218 }
219 }
220 if (dim == 2)
221 {
222 if (xyz_check[0] == dofs.Size())
223 {
224 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 1);
225 }
226 else if (xyz_check[1] == dofs.Size())
227 {
228 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 2);
229 }
230 else
231 {
232 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 4);
233 }
234 }
235 else if (dim == 3)
236 {
237 if (xyz_check[0] == dofs.Size())
238 {
239 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 1);
240 }
241 else if (xyz_check[1] == dofs.Size())
242 {
243 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 2);
244 }
245 else if (xyz_check[2] == dofs.Size())
246 {
247 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 3);
248 }
249 else
250 {
251 pmesh->GetNodalFESpace()->GetMesh()->SetBdrAttribute(i, 4);
252 }
253 }
254 }
255}
256
258 int attr_to_switch)
259{
261 // Switch attribute if all but 1 of the faces of an element will be marked?
262 Array<int> element_attr(pmesh->GetNE());
263 element_attr = 0;
264 for (int e = 0; e < pmesh->GetNE(); e++)
265 {
266 Array<int> faces, ori;
267 if (pmesh->Dimension() == 2)
268 {
269 pmesh->GetElementEdges(e, faces, ori);
270 }
271 else
272 {
273 pmesh->GetElementFaces(e, faces, ori);
274 }
275 int inf1, inf2;
276 int elem1, elem2;
277 int diff_attr_count = 0;
278 int attr1;
279 int attr2;
280 attr1 = mat(e);
281 bool bdr_element = false;
282 element_attr[e] = attr1;
283 int target_attr = -1;
284 for (int f = 0; f < faces.Size(); f++)
285 {
286 pmesh->GetFaceElements(faces[f], &elem1, &elem2);
287 if (elem2 >= 0)
288 {
289 attr2 = elem1 == e ? (int)(mat(elem2)) : (int)(mat(elem1));
290 if (attr1 != attr2 && attr1 == attr_to_switch)
291 {
292 diff_attr_count += 1;
293 target_attr = attr2;
294 }
295 }
296 else
297 {
298 pmesh->GetFaceInfos(faces[f], &inf1, &inf2);
299 if (inf2 >= 0)
300 {
301 Vector dof_vals;
302 Array<int> dofs;
303 mat.GetElementDofValues(pmesh->GetNE() + (-1-elem2), dof_vals);
304 attr2 = (int)(dof_vals(0));
305 if (attr1 != attr2 && attr1 == attr_to_switch)
306 {
307 diff_attr_count += 1;
308 target_attr = attr2;
309 }
310 }
311 else
312 {
313 bdr_element = true;
314 }
315 }
316 }
317
318 if (diff_attr_count == faces.Size()-1 && !bdr_element)
319 {
320 element_attr[e] = target_attr;
321 }
322 }
323 for (int e = 0; e < pmesh->GetNE(); e++)
324 {
325 mat(e) = element_attr[e];
326 pmesh->SetAttribute(e, element_attr[e]+1);
327 }
329 pmesh->SetAttributes();
330}
331
333 FunctionCoefficient &ls_coeff,
334 int amr_iter,
335 ParGridFunction &distance_s,
336 const int quad_order = 5,
337 Array<ParGridFunction *> *pgf_to_update = NULL)
338{
339 mfem::H1_FECollection h1fec(distance_s.ParFESpace()->FEColl()->GetOrder(),
340 pmesh.Dimension());
341 mfem::ParFiniteElementSpace h1fespace(&pmesh, &h1fec);
342 mfem::ParGridFunction x(&h1fespace);
343
344 mfem::L2_FECollection l2fec(0, pmesh.Dimension());
345 mfem::ParFiniteElementSpace l2fespace(&pmesh, &l2fec);
346 mfem::ParGridFunction el_to_refine(&l2fespace);
347
348 mfem::H1_FECollection lhfec(1, pmesh.Dimension());
349 mfem::ParFiniteElementSpace lhfespace(&pmesh, &lhfec);
350 mfem::ParGridFunction lhx(&lhfespace);
351
352 x.ProjectCoefficient(ls_coeff);
354
356 for (int iter = 0; iter < amr_iter; iter++)
357 {
358 el_to_refine = 0.0;
359 for (int e = 0; e < pmesh.GetNE(); e++)
360 {
361 Array<int> dofs;
362 Vector x_vals;
364 h1fespace.GetElementDofs(e, dofs);
365 const IntegrationRule &ir = irRules.Get(pmesh.GetElementGeometry(e),
367 x.GetValues(e, ir, x_vals);
368 real_t min_val = x_vals.Min();
369 real_t max_val = x_vals.Max();
370 // If the zero level set cuts the elements, mark it for refinement
371 if (min_val < 0 && max_val >= 0)
372 {
373 el_to_refine(e) = 1.0;
374 }
375 }
376
377 // Refine an element if its neighbor will be refined
378 for (int inner_iter = 0; inner_iter < 2; inner_iter++)
379 {
380 el_to_refine.ExchangeFaceNbrData();
381 GridFunctionCoefficient field_in_dg(&el_to_refine);
383 for (int e = 0; e < pmesh.GetNE(); e++)
384 {
385 Array<int> dofs;
386 Vector x_vals;
387 lhfespace.GetElementDofs(e, dofs);
388 const IntegrationRule &ir =
390 lhx.GetValues(e, ir, x_vals);
391 real_t max_val = x_vals.Max();
392 if (max_val > 0)
393 {
394 el_to_refine(e) = 1.0;
395 }
396 }
397 }
398
399 // Make the list of elements to be refined
400 Array<int> el_to_refine_list;
401 for (int e = 0; e < el_to_refine.Size(); e++)
402 {
403 if (el_to_refine(e) > 0.0)
404 {
405 el_to_refine_list.Append(e);
406 }
407 }
408
409 int loc_count = el_to_refine_list.Size();
410 int glob_count = loc_count;
411 MPI_Allreduce(&loc_count, &glob_count, 1, MPI_INT, MPI_SUM,
412 pmesh.GetComm());
413 MPI_Barrier(pmesh.GetComm());
414 if (glob_count > 0)
415 {
416 pmesh.GeneralRefinement(el_to_refine_list, 1);
417 }
418
419 // Update
420 h1fespace.Update();
421 x.Update();
422 x.ProjectCoefficient(ls_coeff);
423
424 l2fespace.Update();
425 el_to_refine.Update();
426
427 lhfespace.Update();
428 lhx.Update();
429
430 distance_s.ParFESpace()->Update();
431 distance_s.Update();
432
433 if (pgf_to_update != NULL)
434 {
435 for (int i = 0; i < pgf_to_update->Size(); i++)
436 {
437 (*pgf_to_update)[i]->ParFESpace()->Update();
438 (*pgf_to_update)[i]->Update();
439 }
440 }
441 }
442}
443
445 FunctionCoefficient &ls_coeff,
446 ParGridFunction &distance_s,
447 const int nDiffuse = 2,
448 const int pLapOrder = 5,
449 const int pLapNewton = 50)
450{
451 mfem::H1_FECollection h1fec(distance_s.ParFESpace()->FEColl()->GetOrder(),
452 pmesh.Dimension());
453 mfem::ParFiniteElementSpace h1fespace(&pmesh, &h1fec);
454 mfem::ParGridFunction x(&h1fespace);
455
456 x.ProjectCoefficient(ls_coeff);
458
459 //Now determine distance
460 const real_t dx = AvgElementSize(pmesh);
461 PLapDistanceSolver dist_solver(pLapOrder, pLapNewton);
462
463 ParFiniteElementSpace pfes_s(*distance_s.ParFESpace());
464
465 // Smooth-out Gibbs oscillations from the input level set. The smoothing
466 // parameter here is specified to be mesh dependent with length scale dx.
467 ParGridFunction filt_gf(&pfes_s);
468 PDEFilter filter(pmesh, 1.0 * dx);
469 filter.Filter(ls_coeff, filt_gf);
470 GridFunctionCoefficient ls_filt_coeff(&filt_gf);
471
472 dist_solver.ComputeScalarDistance(ls_filt_coeff, distance_s);
473 distance_s.SetTrueVector();
474 distance_s.SetFromTrueVector();
475
476 DiffuseField(distance_s, nDiffuse);
477 distance_s.SetTrueVector();
478 distance_s.SetFromTrueVector();
479}
480#endif
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition: array.hpp:769
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:727
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:559
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
Definition: fespace.cpp:2950
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:249
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:521
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:150
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:260
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
Container class for integration rules.
Definition: intrules.hpp:416
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:1006
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1439
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:10558
Geometry::Type GetElementGeometry(int i) const
Definition: mesh.hpp:1371
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:6206
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1336
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
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1433
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:7201
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition: mesh.cpp:6985
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
Definition: mesh.hpp:1342
Abstract parallel finite element space.
Definition: pfespace.hpp:29
void GetElementDofs(int i, Array< int > &dofs, DofTransformation &doftrans) const override
The same as GetElementDofs(), but with a user-allocated DofTransformation object. doftrans must be al...
Definition: pfespace.cpp:468
void Update(bool want_transform=true) override
Definition: pfespace.cpp:3414
Class for parallel grid function.
Definition: pgridfunc.hpp:33
void ProjectDiscCoefficient(VectorCoefficient &coeff) override
Project a discontinuous vector coefficient as a grid function on a continuous finite element space....
Definition: pgridfunc.cpp:584
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:562
void Update() override
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
void GetElementDofValues(int el, Vector &dof_vals) const override
Definition: pgridfunc.cpp:543
Class for parallel meshes.
Definition: pmesh.hpp:34
MPI_Comm GetComm() const
Definition: pmesh.hpp:402
void SetAttributes() override
Determine the sets of unique attribute values in domain and boundary elements.
Definition: pmesh.cpp:1589
Vector data type.
Definition: vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:832
real_t Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:922
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
real_t Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1212
void Filter(ParGridFunction &func, ParGridFunction &ffield)
void ComputeScalarDistance(Coefficient &func, ParGridFunction &fdist)
Definition: distance.cpp:107
int dim
Definition: ex24.cpp:53
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
real_t f(const Vector &p)
real_t in_trapezium(const Vector &x, real_t a, real_t b, real_t l)
real_t reactor(const Vector &x)
real_t in_pipe(const Vector &x, int pipedir, Vector x_pipe_center, real_t radius, real_t minv, real_t maxv)
real_t in_cube(const Vector &x, real_t xc, real_t yc, real_t zc, real_t lx, real_t ly, real_t lz)
real_t in_circle(const Vector &x, const Vector &x_center, real_t radius)
real_t in_rectangle(const Vector &x, real_t xc, real_t yc, real_t w, real_t h)
real_t r_intersect(real_t r1, real_t r2)
real_t csg_cubecylsph(const Vector &x)
real_t r_remove(real_t r1, real_t r2)
void ModifyBoundaryAttributesForNodeMovement(ParMesh *pmesh, ParGridFunction &x)
void ModifyAttributeForMarkingDOFS(ParMesh *pmesh, ParGridFunction &mat, int attr_to_switch)
real_t in_parabola(const Vector &x, real_t h, real_t k, real_t t)
void OptimizeMeshWithAMRAroundZeroLevelSet(ParMesh &pmesh, FunctionCoefficient &ls_coeff, int amr_iter, ParGridFunction &distance_s, const int quad_order=5, Array< ParGridFunction * > *pgf_to_update=NULL)
void ComputeScalarDistanceFromLevelSet(ParMesh &pmesh, FunctionCoefficient &ls_coeff, ParGridFunction &distance_s, const int nDiffuse=2, const int pLapOrder=5, const int pLapNewton=50)
real_t circle_level_set(const Vector &x)
real_t r_union(real_t r1, real_t r2)
real_t AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:47
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:22
float real_t
Definition: config.hpp:43
RefCoord t[3]