MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
shaper.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// --------------------------------------------------------------
13// Shaper Miniapp: Resolve material interfaces by mesh refinement
14// --------------------------------------------------------------
15//
16// This miniapp performs multiple levels of adaptive mesh refinement to resolve
17// the interfaces between different "materials" in the mesh, as specified by the
18// given material() function. It can be used as a simple initial mesh generator,
19// for example in the case when the interface is too complex to describe without
20// local refinement. Both conforming and non-conforming refinements are supported.
21//
22// Two additional versions of this miniapp can be found in the miniapps/toys
23// directory: Mandel uses the Shaper algorithm for fractal visualization, while
24// Mondrian convert an image to an AMR mesh suitable for MFEM computations.
25//
26// Compile with: make shaper
27//
28// Sample runs: shaper
29// shaper -m ../../data/inline-tri.mesh
30// shaper -m ../../data/inline-hex.mesh
31// shaper -m ../../data/inline-tet.mesh
32// shaper -m ../../data/amr-quad.mesh
33// shaper -m ../../data/beam-quad.mesh -a -ncl -1 -sd 4
34// shaper -m ../../data/ball-nurbs.mesh
35// shaper -m ../../data/mobius-strip.mesh
36// shaper -m ../../data/square-disc-surf.mesh
37// shaper -m ../../data/star-q3.mesh -sd 2 -ncl -1
38// shaper -m ../../data/fichera-amr.mesh -a -ncl -1
39
40#include "mfem.hpp"
41#include <fstream>
42#include <iostream>
43
44using namespace mfem;
45using namespace std;
46
47// Given a point x, return its material id as an integer. The ids should be
48// positive. If the point is exactly on the interface, return 0.
49//
50// This particular implementation, rescales the mesh to [-1,1]^sdim given the
51// xmin/xmax bounding box, and shapes in a simple annulus/shell with respect to
52// the rescaled coordinates.
53int material(Vector &x, Vector &xmin, Vector &xmax)
54{
55 static real_t p = 2.0;
56
57 // Rescaling to [-1,1]^sdim
58 for (int i = 0; i < x.Size(); i++)
59 {
60 x(i) = (2*x(i)-xmin(i)-xmax(i))/(xmax(i)-xmin(i));
61 }
62
63 // A simple annulus/shell
64 if (x.Normlp(p) > 0.4 && x.Normlp(p) < 0.6) { return 1; }
65 if (x.Normlp(p) < 0.4 || x.Normlp(p) > 0.6) { return 2; }
66 return 0;
67}
68
69int main(int argc, char *argv[])
70{
71 int sd = 2;
72 int nclimit = 1;
73 const char *mesh_file = "../../data/inline-quad.mesh";
74 bool aniso = false;
75
76 // Parse command line
77 OptionsParser args(argc, argv);
78 args.AddOption(&mesh_file, "-m", "--mesh",
79 "Input mesh file to shape materials in.");
80 args.AddOption(&sd, "-sd", "--sub-divisions",
81 "Number of element subdivisions for interface detection.");
82 args.AddOption(&nclimit, "-ncl", "--nc-limit",
83 "Level of hanging nodes allowed (-1 = unlimited).");
84 args.AddOption(&aniso, "-a", "--aniso", "-i", "--iso",
85 "Enable anisotropic refinement of quads and hexes.");
86 args.Parse();
87 if (!args.Good()) { args.PrintUsage(cout); return 1; }
88 args.PrintOptions(cout);
89
90 // Read initial mesh, get dimensions and bounding box
91 Mesh mesh(mesh_file, 1, 1);
92 int dim = mesh.Dimension();
93 int sdim = mesh.SpaceDimension();
94 Vector xmin, xmax;
95 mesh.GetBoundingBox(xmin, xmax);
96
97 // NURBS meshes don't support non-conforming refinement for now
98 if (mesh.NURBSext) { mesh.SetCurvature(2); }
99
100 // Anisotropic refinement not supported for simplex meshes.
101 if (mesh.MeshGenerator() & 1) { aniso = false; }
102
103 // Mesh attributes will be visualized as piece-wise constants
104 L2_FECollection attr_fec(0, dim);
105 FiniteElementSpace attr_fespace(&mesh, &attr_fec);
106 GridFunction attr(&attr_fespace);
107
108 // GLVis server to visualize to
109 char vishost[] = "localhost";
110 int visport = 19916;
111 socketstream sol_sock(vishost, visport);
112 sol_sock.precision(8);
113
114 // Shaping loop
115 for (int iter = 0; 1; iter++)
116 {
118 for (int i = 0; i < mesh.GetNE(); i++)
119 {
120 bool refine = false;
121
122 // Sample materials in each element using "sd" sub-divisions
123 Vector pt;
126 RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, sd, 1);
127 IntegrationRule &ir = RefG->RefPts;
128
129 // Refine any element where different materials are detected. A more
130 // sophisticated logic can be implemented here -- e.g. don't refine
131 // the interfaces between certain materials.
132 Array<int> mat(ir.GetNPoints());
133 real_t matsum = 0.0;
134 for (int j = 0; j < ir.GetNPoints(); j++)
135 {
136 T->Transform(ir.IntPoint(j), pt);
137 int m = material(pt, xmin, xmax);
138 mat[j] = m;
139 matsum += m;
140 if ((int)matsum != m*(j+1))
141 {
142 refine = true;
143 }
144 }
145
146 // Set the element attribute as the "average". Other choices are
147 // possible here too, e.g. attr(i) = mat;
148 attr(i) = round(matsum/ir.GetNPoints());
149
150 // Mark the element for refinement
151 if (refine)
152 {
153 int type = 7;
154 if (aniso)
155 {
156 // Determine the XYZ bitmask for anisotropic refinement.
157 int dx = 0, dy = 0, dz = 0;
158 const int s = sd+1;
159 if (dim == 2)
160 {
161 for (int jj = 0; jj <= sd; jj++)
162 for (int ii = 0; ii < sd; ii++)
163 {
164 dx += abs(mat[jj*s + ii+1] - mat[jj*s + ii]);
165 dy += abs(mat[(ii+1)*s + jj] - mat[ii*s + jj]);
166 }
167 }
168 else if (dim == 3)
169 {
170 for (int kk = 0; kk <= sd; kk++)
171 for (int jj = 0; jj <= sd; jj++)
172 for (int ii = 0; ii < sd; ii++)
173 {
174 dx += abs(mat[(kk*s + jj)*s + ii+1] - mat[(kk*s + jj)*s + ii]);
175 dy += abs(mat[(kk*s + ii+1)*s + jj] - mat[(kk*s + ii)*s + jj]);
176 dz += abs(mat[((ii+1)*s + jj)*s + kk] - mat[(ii*s + jj)*s + kk]);
177 }
178 }
179 type = 0;
180 const int tol = mat.Size() / 10;
181 if (dx > tol) { type |= 1; }
182 if (dy > tol) { type |= 2; }
183 if (dz > tol) { type |= 4; }
184 if (!type) { type = 7; } // because of tol
185 }
186
187 refs.Append(Refinement(i, type));
188 }
189 }
190
191 // Visualization
192 sol_sock << "solution\n" << mesh << attr;
193 if (iter == 0 && sdim == 2)
194 {
195 sol_sock << "keys 'RjlmpppppppppppppA*************'\n";
196 }
197 if (iter == 0 && sdim == 3)
198 {
199 sol_sock << "keys 'YYYYYYYYYXXXXXXXmA********8888888pppttt";
200 if (dim == 3) { sol_sock << "iiM"; }
201 sol_sock << "'\n";
202 }
203 sol_sock << flush;
204
205 // Ask the user if we should continue refining
206 char yn;
207 cout << "Mesh has " << mesh.GetNE() << " elements. \n"
208 << "Continue shaping? --> ";
209 cin >> yn;
210 if (yn == 'n' || yn == 'q') { break; }
211
212 // Perform refinement, update spaces and grid functions
213 mesh.GeneralRefinement(refs, -1, nclimit);
214 attr_fespace.Update();
215 attr.Update();
216 }
217
218 // Set element attributes in the mesh object before saving
219 for (int i = 0; i < mesh.GetNE(); i++)
220 {
221 mesh.SetAttribute(i, static_cast<int>(attr(i)));
222 }
223 mesh.SetAttributes();
224
225 // Save the final mesh
226 ofstream mesh_ofs("shaper.mesh");
227 mesh_ofs.precision(8);
228 mesh.Print(mesh_ofs);
229}
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
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3439
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1136
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:259
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Mesh data type.
Definition: mesh.hpp:56
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:290
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:10558
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1336
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition: mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:138
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition: mesh.cpp:357
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1163
int MeshGenerator() const
Get the mesh generator/type.
Definition: mesh.hpp:1187
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:6211
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1385
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition: mesh.cpp:1604
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
IntegrationRule RefPts
Definition: geom.hpp:317
Vector data type.
Definition: vector.hpp:80
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
real_t Normlp(real_t p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:872
int dim
Definition: ex24.cpp:53
int visport
char vishost[]
int main()
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1891
float real_t
Definition: config.hpp:43
real_t p(const Vector &x, real_t t)
Definition: navier_mms.cpp:53
RefCoord s[3]
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition: shaper.cpp:53