MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
fit-node-position.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// Fitting of Selected Mesh Nodes to Specified Physical Positions
14// ------------------------------------------------------------------
15//
16// This example fits a selected set of the mesh nodes to given physical
17// positions while maintaining a valid mesh with good quality.
18//
19// Sample runs:
20// mpirun -np 4 fit-node-position
21// mpirun -np 4 fit-node-position -m square01-tri.mesh
22// mpirun -np 4 fit-node-position -m ./cube.mesh
23// mpirun -np 4 fit-node-position -m ./cube-tet.mesh -rs 0
24
25#include "mfem.hpp"
26#include "../common/mfem-common.hpp"
27
28using namespace mfem;
29using namespace std;
30
31char vishost[] = "localhost";
32int visport = 19916;
33int wsize = 350;
34
35int main (int argc, char *argv[])
36{
37 // Initialize MPI.
38 Mpi::Init();
39 int myid = Mpi::WorldRank();
40
41 const char *mesh_file = "square01.mesh";
42 int rs_levels = 2;
43 int mesh_poly_deg = 2;
44 int quad_order = 5;
45 bool glvis = true;
46
47 // Parse command-line options.
48 OptionsParser args(argc, argv);
49 args.AddOption(&mesh_file, "-m", "--mesh",
50 "Mesh file to use.");
51 args.AddOption(&rs_levels, "-rs", "--refine-serial",
52 "Number of times to refine the mesh uniformly in serial.");
53 args.AddOption(&mesh_poly_deg, "-o", "--order",
54 "Polynomial degree of mesh finite element space.");
55 args.AddOption(&quad_order, "-qo", "--quad_order",
56 "Order of the quadrature rule.");
57 args.AddOption(&glvis, "-vis", "--visualization", "-no-vis",
58 "--no-visualization",
59 "Enable or disable GLVis visualization.");
60 args.Parse();
61 if (!args.Good())
62 {
63 if (myid == 0) { args.PrintUsage(cout); }
64 return 1;
65 }
66 if (myid == 0) { args.PrintOptions(cout); }
67
68 // Read and refine the mesh.
69 Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
70 for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
71 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
72 delete mesh;
73 const int dim = pmesh.Dimension();
74
75 // Setup mesh curvature and GridFunction that stores the coordinates.
77 if (mesh_poly_deg <= 0)
78 {
79 fec_mesh = new QuadraticPosFECollection;
80 mesh_poly_deg = 2;
81 }
82 else { fec_mesh = new H1_FECollection(mesh_poly_deg, dim); }
83 ParFiniteElementSpace pfes_mesh(&pmesh, fec_mesh, dim);
84 pmesh.SetNodalFESpace(&pfes_mesh);
85 ParGridFunction coord(&pfes_mesh);
86 pmesh.SetNodalGridFunction(&coord);
87 ParGridFunction x0(coord);
88
89 // Pick which nodes to fit and select the target positions.
90 // (attribute 2 would have a prescribed deformation in y-direction, same x).
91 Array<bool> fit_marker(pfes_mesh.GetNDofs());
92 ParGridFunction fit_marker_vis_gf(&pfes_mesh);
93 ParGridFunction coord_target(&pfes_mesh);
94 Array<int> vdofs;
95 fit_marker = false;
96 coord_target = coord;
97 fit_marker_vis_gf = 0.0;
98 for (int e = 0; e < pmesh.GetNBE(); e++)
99 {
100 const int nd = pfes_mesh.GetBE(e)->GetDof();
101 const int attr = pmesh.GetBdrElement(e)->GetAttribute();
102 if (attr != 2) { continue; }
103
104 pfes_mesh.GetBdrElementVDofs(e, vdofs);
105 for (int j = 0; j < nd; j++)
106 {
107 int j_x = vdofs[j], j_y = vdofs[nd+j];
108 const real_t x = coord(j_x),
109 z = (dim == 2) ? 0.0 : coord(vdofs[2*nd + j]);
110 fit_marker[pfes_mesh.VDofToDof(j_x)] = true;
111 fit_marker_vis_gf(j_x) = 1.0;
112 if (coord(j_y) < 0.5)
113 {
114 coord_target(j_y) = 0.1 * sin(4 * M_PI * x) * cos(M_PI * z);
115 }
116 else
117 {
118 if (coord(j_x) < 0.5)
119 {
120 coord_target(j_y) = 1.0 + 0.1 * sin(2 * M_PI * x);
121 }
122 else
123 {
124 coord_target(j_y) = 1.0 + 0.1 * sin(2 * M_PI * (x + 0.5));
125 }
126
127 }
128 }
129 }
130
131 // Visualize the selected nodes and their target positions.
132 if (glvis)
133 {
134 socketstream vis1;
135 coord = coord_target;
136 common::VisualizeField(vis1, "localhost", 19916, fit_marker_vis_gf,
137 "Target positions (DOFS with value 1)",
138 0, 0, 400, 400, (dim == 2) ? "Rjm" : "");
139 coord = x0;
140 }
141
142 // Allow slipping along the remaining boundaries.
143 // (attributes 1 and 3 would slip, while 4 is completely fixed).
144 int n = 0;
145 for (int i = 0; i < pmesh.GetNBE(); i++)
146 {
147 const int nd = pfes_mesh.GetBE(i)->GetDof();
148 const int attr = pmesh.GetBdrElement(i)->GetAttribute();
149 MFEM_VERIFY(!(dim == 2 && attr == 3),
150 "Boundary attribute 3 must be used only for 3D meshes. "
151 "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
152 "components, rest for free nodes), or use -fix-bnd.");
153 if (attr == 1 || attr == 3) { n += nd; }
154 if (attr == 4) { n += nd * dim; }
155 }
156 Array<int> ess_vdofs(n);
157 n = 0;
158 for (int i = 0; i < pmesh.GetNBE(); i++)
159 {
160 const int nd = pfes_mesh.GetBE(i)->GetDof();
161 const int attr = pmesh.GetBdrElement(i)->GetAttribute();
162 pfes_mesh.GetBdrElementVDofs(i, vdofs);
163 if (attr == 1) // Fix x components.
164 {
165 for (int j = 0; j < nd; j++)
166 { ess_vdofs[n++] = vdofs[j]; }
167 }
168 else if (attr == 3) // Fix z components.
169 {
170 for (int j = 0; j < nd; j++)
171 { ess_vdofs[n++] = vdofs[j+2*nd]; }
172 }
173 else if (attr == 4) // Fix all components.
174 {
175 for (int j = 0; j < vdofs.Size(); j++)
176 { ess_vdofs[n++] = vdofs[j]; }
177 }
178 }
179
180 // TMOP setup.
181 TMOP_QualityMetric *metric;
182 if (dim == 2) { metric = new TMOP_Metric_002; }
183 else { metric = new TMOP_Metric_302; }
185 pfes_mesh.GetComm());
186 ConstantCoefficient fit_weight(100.0);
187 auto integ = new TMOP_Integrator(metric, &target, nullptr);
188 integ->EnableSurfaceFitting(coord_target, fit_marker, fit_weight);
189
190 // Linear solver.
191 MINRESSolver minres(pfes_mesh.GetComm());
192 minres.SetMaxIter(100);
193 minres.SetRelTol(1e-12);
194 minres.SetAbsTol(0.0);
195
196 // Nonlinear solver.
197 ParNonlinearForm a(&pfes_mesh);
198 a.SetEssentialVDofs(ess_vdofs);
199 a.AddDomainIntegrator(integ);
200 const IntegrationRule &ir =
201 IntRules.Get(pfes_mesh.GetFE(0)->GetGeomType(), quad_order);
202 TMOPNewtonSolver solver(pfes_mesh.GetComm(), ir, 0);
203 solver.SetOperator(a);
204 solver.SetPreconditioner(minres);
205 solver.SetPrintLevel(1);
206 solver.SetMaxIter(200);
207 solver.SetRelTol(1e-10);
208 solver.SetAbsTol(0.0);
211
212 // Solve.
213 Vector b(0);
214 coord.SetTrueVector();
215 solver.Mult(b, coord.GetTrueVector());
216 coord.SetFromTrueVector();
217
218 if (glvis)
219 {
220 socketstream vis2;
221 common::VisualizeMesh(vis2, "localhost", 19916, pmesh, "Final mesh",
222 400, 0, 400, 400);
223 }
224
225 delete metric;
226 return 0;
227}
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
A coefficient that is constant across space and time.
Definition: coefficient.hpp:85
int GetAttribute() const
Return element's attribute.
Definition: element.hpp:55
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:27
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition: fespace.cpp:3204
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:710
int VDofToDof(int vdof) const
Compute the inverse of the Dof to VDof mapping for a single index vdof.
Definition: fespace.hpp:995
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
Definition: fespace.cpp:303
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:326
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:329
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
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:130
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
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:1006
void SetRelTol(real_t rtol)
Definition: solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void SetMaxIter(int max_it)
Definition: solvers.hpp:211
void SetAbsTol(real_t atol)
Definition: solvers.hpp:210
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:787
MINRES method.
Definition: solvers.hpp:628
Mesh data type.
Definition: mesh.hpp:56
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition: mesh.hpp:1298
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:6200
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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
Abstract parallel finite element space.
Definition: pfespace.hpp:29
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
const FiniteElement * GetFE(int i) const override
Definition: pfespace.cpp:534
Class for parallel grid function.
Definition: pgridfunc.hpp:33
Class for parallel meshes.
Definition: pmesh.hpp:34
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2028
Parallel non-linear operator on the true dofs.
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:796
void EnableAdaptiveSurfaceFitting()
Definition: tmop_tools.hpp:235
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:261
void SetTerminationWithMaxSurfaceFittingError(real_t max_error)
Definition: tmop_tools.hpp:252
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:286
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1738
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:687
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P....
Definition: tmop.hpp:24
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1334
Vector data type.
Definition: vector.hpp:80
int dim
Definition: ex24.cpp:53
int wsize
int visport
char vishost[]
int main()
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
void VisualizeMesh(socketstream &sock, const char *vishost, int visport, Mesh &mesh, const char *title, int x, int y, int w, int h, const char *keys)
Definition: fem_extras.cpp:59
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
Definition: fem_extras.cpp:92
float real_t
Definition: config.hpp:43
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:486