MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex1p.cpp
Go to the documentation of this file.
1 // MFEM + Moonolith Example 1 (parallel version)
2 //
3 // Compile with: make ex1p
4 //
5 // Moonolith sample runs:
6 // mpirun -np 4 ex1p
7 // mpirun -np 4 ex1p --source_refinements 1 --dest_refinements 2
8 // mpirun -np 4 ex1p -s ../../data/inline-hex.mesh -d ../../data/inline-tet.mesh
9 //
10 // Description: This example code demonstrates the use of MFEM for transferring
11 // discrete fields from one finite element mesh to another. The
12 // meshes can be of arbitrary shape and completely unrelated with
13 // each other. This feature can be used for implementing immersed
14 // domain methods for fluid-structure interaction or general
15 // multi-physics applications.
16 //
17 // This particular example is for parallel runtimes. Vector FE is
18 // an experimental feature in parallel.
19 
20 #include "example_utils.hpp"
21 #include "mfem.hpp"
22 
23 using namespace mfem;
24 using namespace std;
25 
26 void destination_transform(const Vector &x, Vector &x_new)
27 {
28  x_new = x;
29  // x_new *= 0.5;
30 }
31 
32 int main(int argc, char *argv[])
33 {
34  MPI_Init(&argc, &argv);
35 
36  int num_procs, rank;
37 
38  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
39  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
40 
41  // Init transfer library context, with MPI handled outside the library
42  InitTransfer(argc, argv, MPI_COMM_WORLD);
43 
44  const char *source_mesh_file = "../../data/inline-tri.mesh";
45  const char *destination_mesh_file = "../../data/inline-quad.mesh";
46 
47  int src_n_refinements = 0;
48  int dest_n_refinements = 0;
49  int source_fe_order = 1;
50  int dest_fe_order = 1;
51  bool visualization = true;
52  bool use_vector_fe = false;
53  bool verbose = false;
54  bool assemble_mass_and_coupling_together = true;
55 
56  OptionsParser args(argc, argv);
57  args.AddOption(&source_mesh_file, "-s", "--source_mesh",
58  "Mesh file to use for src.");
59  args.AddOption(&destination_mesh_file, "-d", "--destination_mesh",
60  "Mesh file to use for dest.");
61  args.AddOption(&src_n_refinements, "-sr", "--source_refinements",
62  "Number of src refinements");
63  args.AddOption(&dest_n_refinements, "-dr", "--dest_refinements",
64  "Number of dest refinements");
65  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
66  "--no-visualization",
67  "Enable or disable GLVis visualization.");
68  args.AddOption(&source_fe_order, "-so", "--source_fe_order",
69  "Order of the src finite elements");
70  args.AddOption(&dest_fe_order, "-do", "--dest_fe_order",
71  "Order of the dest finite elements");
72  args.AddOption(&verbose, "-verb", "--verbose", "--no-verb", "--no-verbose",
73  "Enable/Disable verbose output");
74  args.AddOption(&use_vector_fe, "-vfe", "--use_vector_fe", "-no-vfe",
75  "--no-vector_fe", "Use vector finite elements (Experimental)");
76  args.AddOption(&assemble_mass_and_coupling_together, "-act",
77  "--assemble_mass_and_coupling_together", "-no-act",
78  "--no-assemble_mass_and_coupling_together",
79  "Assemble mass and coupling operators together (better for non-affine elements)");
80  args.Parse();
81  check_options(args);
82 
83  shared_ptr<Mesh> src_mesh, dest_mesh;
84 
85  ifstream imesh;
86 
87  imesh.open(destination_mesh_file);
88  if (imesh)
89  {
90  dest_mesh = make_shared<Mesh>(imesh, 1, 1);
91  imesh.close();
92  }
93  else
94  {
95  if (rank == 0)
96  mfem::err << "WARNING: Destination mesh file not found: "
97  << destination_mesh_file << "\n"
98  << "Using default 2D quad mesh.";
99 
100  dest_mesh = make_shared<Mesh>(4, 4, Element::QUADRILATERAL);
101  }
102 
103  const int dim = dest_mesh->Dimension();
104 
105  dest_mesh->Transform(&destination_transform);
106 
107  Vector box_min(dim), box_max(dim), range(dim);
108  dest_mesh->GetBoundingBox(box_min, box_max);
109  range = box_max;
110  range -= box_min;
111 
112  imesh.open(source_mesh_file);
113 
114  if (imesh)
115  {
116  src_mesh = make_shared<Mesh>(imesh, 1, 1);
117  imesh.close();
118  }
119  else
120  {
121  if (rank == 0)
122  mfem::err << "WARNING: Source mesh file not found: " << source_mesh_file
123  << "\n"
124  << "Using default box mesh.\n";
125 
126  if (dim == 2)
127  {
128  src_mesh =
129  make_shared<Mesh>(4, 4, Element::TRIANGLE, 1, range[0], range[1]);
130  }
131  else if (dim == 3)
132  {
133  src_mesh = make_shared<Mesh>(4, 4, 4, Element::TETRAHEDRON, 1, range[0],
134  range[1], range[2]);
135  }
136 
137  for (int i = 0; i < src_mesh->GetNV(); ++i)
138  {
139  double *v = src_mesh->GetVertex(i);
140 
141  for (int d = 0; d < dim; ++d)
142  {
143  v[d] += box_min[d];
144  }
145  }
146  }
147 
148  for (int i = 0; i < src_n_refinements; ++i)
149  {
150  src_mesh->UniformRefinement();
151  }
152 
153  for (int i = 0; i < dest_n_refinements; ++i)
154  {
155  dest_mesh->UniformRefinement();
156  }
157 
158  auto p_src_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *src_mesh);
159  auto p_dest_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *dest_mesh);
160 
161  shared_ptr<FiniteElementCollection> src_fe_coll, dest_fe_coll;
162 
163  if (use_vector_fe)
164  {
165  src_fe_coll =
166  make_shared<RT_FECollection>(source_fe_order, src_mesh->Dimension());
167  dest_fe_coll =
168  make_shared<RT_FECollection>(dest_fe_order, dest_mesh->Dimension());
169  }
170  else
171  {
172  src_fe_coll =
173  make_shared<L2_FECollection>(source_fe_order, src_mesh->Dimension());
174  dest_fe_coll =
175  make_shared<L2_FECollection>(dest_fe_order, dest_mesh->Dimension());
176  }
177 
178  auto src_fe =
179  make_shared<ParFiniteElementSpace>(p_src_mesh.get(), src_fe_coll.get());
180 
181  auto dest_fe =
182  make_shared<ParFiniteElementSpace>(p_dest_mesh.get(), dest_fe_coll.get());
183 
184  ParGridFunction src_fun(src_fe.get());
185 
186  // To be used with standard fe
188 
189  // To be used with vector fe
190  VectorFunctionCoefficient vector_coeff(dim, &vector_fun);
191 
192  if (use_vector_fe)
193  {
194  src_fun.ProjectCoefficient(vector_coeff);
195  src_fun.Update();
196  }
197  else
198  {
199  src_fun.ProjectCoefficient(coeff);
200  src_fun.Update();
201  }
202 
203  ParGridFunction dest_fun(dest_fe.get());
204  dest_fun = 0.0;
205  dest_fun.Update();
206 
207  ParMortarAssembler assembler(src_fe, dest_fe);
209  assemble_mass_and_coupling_together);
210  assembler.SetVerbose(verbose);
211 
212  if (use_vector_fe)
213  {
214  assembler.AddMortarIntegrator(make_shared<VectorL2MortarIntegrator>());
215  }
216  else
217  {
218  assembler.AddMortarIntegrator(make_shared<L2MortarIntegrator>());
219  }
220 
221  if (assembler.Transfer(src_fun, dest_fun))
222  {
223 
224  if (visualization)
225  {
226  double src_err = 0;
227  double dest_err = 0;
228 
229  if (use_vector_fe)
230  {
231  src_err = src_fun.ComputeL2Error(vector_coeff);
232  dest_err = dest_fun.ComputeL2Error(vector_coeff);
233  }
234  else
235  {
236  src_err = src_fun.ComputeL2Error(coeff);
237  dest_err = dest_fun.ComputeL2Error(coeff);
238  }
239 
240  if (rank == 0)
241  {
242  mfem::out << "l2 error: src: " << src_err << ", dest: " << dest_err
243  << std::endl;
244  }
245 
246  plot(*p_src_mesh, src_fun, "source");
247  plot(*p_dest_mesh, dest_fun, "destination");
248  }
249  }
250  else
251  {
252  mfem::out << "No intersection no transfer!" << std::endl;
253  }
254 
255  // Finalize transfer library context
257  return MPI_Finalize();
258 }
void SetVerbose(const bool verbose)
Expose process details with verbose output.
void destination_transform(const Vector &x, Vector &x_new)
Definition: ex1p.cpp:26
This class implements the parallel 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 particular code is also used with LLNL for large scale multilevel Monte Carlo simulations. This algorithm allows to perform quadrature in the intersection of elements of two separate, unrelated, and arbitrarily distributed 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. Convex non-affine elements are partially supported, however, high order (&gt;3) finite element discretizations or nonlinear geometric transformations might generate undesired oscillations. Discontinuous fields in general can only be mapped to order 0 destination fields. For such cases localized versions of the projection will have to be developed.
bool Transfer(const ParGridFunction &src_fun, ParGridFunction &dest_fun)
transfer a function from source to destination. if the transfer is to be performed multiple times use...
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
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...
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
A general vector function coefficient.
void plot(mfem::Mesh &mesh, mfem::GridFunction &x, std::string title)
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 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
A general function coefficient.
void SetAssembleMassAndCouplingTogether(const bool value)
Control if the Mass matrix is computed together with the coupling operator every time.
Vector data type.
Definition: vector.hpp:60
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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
int FinalizeTransfer()
Finalize the par_moonolith library.
Definition: transfer.cpp:23
int main()