MFEM  v4.5.2
Finite element discretization library
ex1.cpp
Go to the documentation of this file.
1 // MFEM + Moonolith Example 1
2 //
3 // Compile with: make ex1
4 //
5 // Moonolith sample runs:
6 // ex1
7 // ex1 --source_refinements 1 --dest_refinements 2
8 // ex1 --source_refinements 1 --dest_refinements 2 --use_vector_fe
9 // ex1 -s ../../data/inline-hex.mesh -d ../../data/inline-tet.mesh
10 //
11 // Description: This example code demonstrates the use of MFEM for transferring
12 // discrete fields from one finite element mesh to another. The
13 // meshes can be of arbitrary shape and completely unrelated with
14 // each other. This feature can be used for implementing immersed
15 // domain methods for fluid-structure interaction or general
16 // multi-physics applications.
17 //
18 // This particular example is only for serial runtimes.
19 
20 #include "example_utils.hpp"
21 #include "mfem.hpp"
22 
23 using namespace mfem;
24 using namespace std;
25 
26 int main(int argc, char *argv[])
27 {
28  // Init transfer library context
29  InitTransfer(argc, argv);
30 
31  const char *source_mesh_file = "../../data/inline-tri.mesh";
32  const char *destination_mesh_file = "../../data/inline-quad.mesh";
33 
34  int src_n_refinements = 0;
35  int dest_n_refinements = 0;
36  int source_fe_order = 1;
37  int dest_fe_order = 1;
38  bool visualization = true;
39  bool use_vector_fe = false;
40  bool verbose = false;
41 
42  OptionsParser args(argc, argv);
43  args.AddOption(&source_mesh_file, "-s", "--source_mesh",
44  "Mesh file to use for src.");
45  args.AddOption(&destination_mesh_file, "-d", "--destination_mesh",
46  "Mesh file to use for dest.");
47  args.AddOption(&src_n_refinements, "-sr", "--source_refinements",
48  "Number of src refinements");
49  args.AddOption(&dest_n_refinements, "-dr", "--dest_refinements",
50  "Number of dest refinements");
51  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
52  "--no-visualization",
53  "Enable or disable GLVis visualization.");
54  args.AddOption(&source_fe_order, "-so", "--source_fe_order",
55  "Order of the src finite elements");
56  args.AddOption(&dest_fe_order, "-do", "--dest_fe_order",
57  "Order of the dest finite elements");
58  args.AddOption(&verbose, "-verb", "--verbose", "--no-verb", "--no-verbose",
59  "Enable/Disable verbose output");
60  args.AddOption(&use_vector_fe, "-vfe", "--use_vector_fe", "-no-vfe",
61  "--no-vector_fe", "Use vector finite elements");
62  args.Parse();
63  check_options(args);
64 
65  shared_ptr<Mesh> src_mesh, dest_mesh;
66 
67  ifstream imesh;
68 
69  imesh.open(destination_mesh_file);
70  if (imesh)
71  {
72  dest_mesh = make_shared<Mesh>(imesh, 1, 1);
73  imesh.close();
74  }
75  else
76  {
77  mfem::err << "WARNING: Destination mesh file not found: "
78  << destination_mesh_file << "\n"
79  << "Using default 2D quad mesh.";
80 
81  dest_mesh = make_shared<Mesh>(4, 4, Element::QUADRILATERAL);
82  }
83 
84  const int dim = dest_mesh->Dimension();
85 
86  Vector box_min(dim), box_max(dim), range(dim);
87  dest_mesh->GetBoundingBox(box_min, box_max);
88  range = box_max;
89  range -= box_min;
90 
91  imesh.open(source_mesh_file);
92 
93  if (imesh)
94  {
95  src_mesh = make_shared<Mesh>(imesh, 1, 1);
96  imesh.close();
97  }
98  else
99  {
100  mfem::err << "WARNING: Source mesh file not found: " << source_mesh_file
101  << "\n"
102  << "Using default box mesh.\n";
103 
104  if (dim == 2)
105  {
106  src_mesh =
107  make_shared<Mesh>(4, 4, Element::TRIANGLE, 1, range[0], range[1]);
108  }
109  else if (dim == 3)
110  {
111  src_mesh = make_shared<Mesh>(4, 4, 4, Element::TETRAHEDRON, 1, range[0],
112  range[1], range[2]);
113  }
114 
115  for (int i = 0; i < src_mesh->GetNV(); ++i)
116  {
117  double *v = src_mesh->GetVertex(i);
118 
119  for (int d = 0; d < dim; ++d)
120  {
121  v[d] += box_min[d];
122  }
123  }
124  }
125 
126  for (int i = 0; i < src_n_refinements; ++i)
127  {
128  src_mesh->UniformRefinement();
129  }
130 
131  for (int i = 0; i < dest_n_refinements; ++i)
132  {
133  dest_mesh->UniformRefinement();
134  }
135 
136  shared_ptr<FiniteElementCollection> src_fe_coll, dest_fe_coll;
137 
138  if (use_vector_fe)
139  {
140  src_fe_coll =
141  make_shared<RT_FECollection>(source_fe_order, src_mesh->Dimension());
142  dest_fe_coll =
143  make_shared<RT_FECollection>(dest_fe_order, dest_mesh->Dimension());
144  }
145  else
146  {
147  src_fe_coll =
148  make_shared<L2_FECollection>(source_fe_order, src_mesh->Dimension());
149  dest_fe_coll =
150  make_shared<L2_FECollection>(dest_fe_order, dest_mesh->Dimension());
151  }
152 
153  auto src_fe =
154  make_shared<FiniteElementSpace>(src_mesh.get(), src_fe_coll.get());
155 
156  auto dest_fe =
157  make_shared<FiniteElementSpace>(dest_mesh.get(), dest_fe_coll.get());
158 
159  GridFunction src_fun(src_fe.get());
160  GridFunction dest_fun(dest_fe.get());
161  src_fun = 1.0;
162 
163  // To be used with standard fe
165 
166  // To be used with vector fe
167  VectorFunctionCoefficient vector_coeff(dim, &vector_fun);
168 
169  if (use_vector_fe)
170  {
171  src_fun.ProjectCoefficient(vector_coeff);
172  src_fun.Update();
173  }
174  else
175  {
176  src_fun.ProjectCoefficient(coeff);
177  src_fun.Update();
178  }
179 
180  dest_fun = 0.0;
181  dest_fun.Update();
182 
183  MortarAssembler assembler(src_fe, dest_fe);
184  assembler.SetVerbose(verbose);
185 
186  if (use_vector_fe)
187  {
188  assembler.AddMortarIntegrator(make_shared<VectorL2MortarIntegrator>());
189  }
190  else
191  {
192  assembler.AddMortarIntegrator(make_shared<L2MortarIntegrator>());
193  }
194 
195  if (assembler.Transfer(src_fun, dest_fun))
196  {
197  if (visualization)
198  {
199  dest_fun.Update();
200 
201  double src_err = 0;
202  double dest_err = 0;
203 
204  if (use_vector_fe)
205  {
206  src_err = src_fun.ComputeL2Error(vector_coeff);
207  dest_err = dest_fun.ComputeL2Error(vector_coeff);
208  }
209  else
210  {
211  src_err = src_fun.ComputeL2Error(coeff);
212  dest_err = dest_fun.ComputeL2Error(coeff);
213  }
214 
215  mfem::out << "l2 error: src: " << src_err << ", dest: " << dest_err
216  << std::endl;
217 
218  plot(*src_mesh, src_fun, "source");
219  plot(*dest_mesh, dest_fun, "destination");
220  }
221  }
222  else
223  {
224  mfem::out << "No intersection -> no transfer!" << std::endl;
225  }
226 
227  // Finalize transfer library context
228  return FinalizeTransfer();
229 }
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void SetVerbose(const bool verbose)
Expose process details with verbose output.
int Dimension() const
Definition: mesh.hpp:1047
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
STL namespace.
This class implements the serial variational transfer between finite element spaces. Variational transfer has been shown to have better approximation properties than standard interpolation. This facilities can be used for supporting applications which require the handling of non matching meshes. For instance: General multi-physics problems, fluid structure interaction, or even visualization of average quantities within subvolumes This algorithm allows to perform quadrature in the intersection of elements of two separate and unrelated meshes. It generates quadrature rules in the intersection which allows us to integrate-with to machine precision using the mfem::MortarIntegrator interface. See https://doi.org/10.1137/15M1008361 for and in-depth explanation. At this time curved elements are not supported.
int main(int argc, char *argv[])
Definition: ex1.cpp:74
void check_options(mfem::OptionsParser &args)
void InitTransfer(int argc, char *argv[])
Initializes the par_moonolith library. It also calls MPI_Init.
Definition: transfer.cpp:19
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:933
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
bool Transfer(const GridFunction &src_fun, GridFunction &dest_fun)
transfer a function from source to destination. if the transfer is to be performed multiple times use...
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1053
A general vector function coefficient.
void plot(mfem::Mesh &mesh, mfem::GridFunction &x, std::string title)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
double example_fun(const mfem::Vector &x)
void vector_fun(const mfem::Vector &x, mfem::Vector &f)
int dim
Definition: ex24.cpp:53
void AddMortarIntegrator(const std::shared_ptr< MortarIntegrator > &integrator)
This method must be called before Assemble or Transfer. It will assemble the operator in all intersec...
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
int FinalizeTransfer()
Finalize the par_moonolith library.
Definition: transfer.cpp:23