MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex2p.cpp
Go to the documentation of this file.
1 // MFEM + Moonolith Example 2 (parallel version)
2 //
3 // Compile with: make ex2p
4 //
5 // Moonolith sample runs:
6 // mpirun -np 4 ex2p
7 // mpirun -np 4 ex2p --source_refinements 1 --dest_refinements 2
8 // mpirun -np 4 ex2p -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 concerns discontinuous Galerkin FEM with
18 // adaptive mesh refinement for parallel runtimes.
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 *= .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 
50  // Source fe order has to be greater or equal than destination order
51  int source_fe_order = 1;
52  int dest_fe_order = 0;
53  bool visualization = true;
54  bool verbose = false;
55  int max_iterations = 30000;
56 
57  OptionsParser args(argc, argv);
58  args.AddOption(&source_mesh_file, "-s", "--source_mesh",
59  "Mesh file to use for src.");
60  args.AddOption(&destination_mesh_file, "-d", "--destination_mesh",
61  "Mesh file to use for dest.");
62  args.AddOption(&src_n_refinements, "-sr", "--source_refinements",
63  "Number of src refinements");
64  args.AddOption(&dest_n_refinements, "-dr", "--dest_refinements",
65  "Number of dest refinements");
66  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
67  "--no-visualization",
68  "Enable or disable GLVis visualization.");
69  args.AddOption(&source_fe_order, "-so", "--source_fe_order",
70  "Order of the src finite elements");
71  args.AddOption(&dest_fe_order, "-do", "--dest_fe_order",
72  "Order of the dest finite elements");
73  args.AddOption(&verbose, "-verb", "--verbose", "--no-verb", "--no-verbose",
74  "Enable/Disable verbose output");
75  args.AddOption(&max_iterations, "-m", "--max_iterations",
76  "Max number of solver iterations");
77  args.Parse();
78  check_options(args);
79 
80  if (source_fe_order == 0 && dest_fe_order != 0)
81  {
82  mfem::out <<
83  "Source fe order should not be 0 unless destination fe order is also 0!\n";
84 
86  return MPI_Finalize();
87  }
88 
89  ifstream imesh(source_mesh_file);
90  shared_ptr<Mesh> src_mesh, dest_mesh;
91  if (imesh)
92  {
93  src_mesh = make_shared<Mesh>(imesh, 1, 1);
94  imesh.close();
95  }
96  else
97  {
98  if (rank == 0)
99  mfem::err << "WARNING: Source mesh file not found: " << source_mesh_file
100  << "\n"
101  << "Using default 2D triangle mesh.";
102  src_mesh = make_shared<Mesh>(4, 4, Element::TRIANGLE);
103  }
104 
105  imesh.open(destination_mesh_file);
106  if (imesh)
107  {
108  dest_mesh = make_shared<Mesh>(imesh, 1, 1);
109  imesh.close();
110  }
111  else
112  {
113  if (rank == 0)
114  mfem::err << "WARNING: Destination mesh file not found: "
115  << destination_mesh_file << "\n"
116  << "Using default 2D quad mesh.";
117  dest_mesh = make_shared<Mesh>(4, 4, Element::QUADRILATERAL);
118  }
119 
120  dest_mesh->Transform(&destination_transform);
121 
122  for (int i = 0; i < src_n_refinements; ++i)
123  {
124  src_mesh->UniformRefinement();
125  }
126 
127  for (int i = 0; i < dest_n_refinements; ++i)
128  {
129  dest_mesh->UniformRefinement();
130  }
131 
132  src_mesh->EnsureNCMesh();
133  dest_mesh->EnsureNCMesh();
134 
135  {
136  for (int l = 0; l < 4; l++)
137  {
138  src_mesh->RandomRefinement(0.1); // 10% probability
139  }
140  }
141 
142  {
143  for (int l = 0; l < 4; l++)
144  {
145  dest_mesh->RandomRefinement(0.1); // 10% probability
146  }
147  }
148 
149  auto p_src_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *src_mesh);
150  auto p_dest_mesh = make_shared<ParMesh>(MPI_COMM_WORLD, *dest_mesh);
151 
152  auto src_fe_coll =
153  make_shared<DG_FECollection>(source_fe_order, p_src_mesh->Dimension());
154  auto src_fe =
155  make_shared<ParFiniteElementSpace>(p_src_mesh.get(), src_fe_coll.get());
156 
157  auto dest_fe_coll =
158  make_shared<DG_FECollection>(dest_fe_order, p_dest_mesh->Dimension());
159  auto dest_fe =
160  make_shared<ParFiniteElementSpace>(p_dest_mesh.get(), dest_fe_coll.get());
161 
162  ParGridFunction src_fun(src_fe.get());
164  make_fun(*src_fe, coeff, src_fun);
165 
166  ParGridFunction dest_fun(dest_fe.get());
167  dest_fun = 0.0;
168  dest_fun.Update();
169 
170  ParMortarAssembler assembler(src_fe, dest_fe);
171  assembler.SetVerbose(verbose);
172  assembler.SetMaxSolverIterations(max_iterations);
173 
174  assembler.AddMortarIntegrator(make_shared<L2MortarIntegrator>());
175  if (assembler.Transfer(src_fun, dest_fun))
176  {
177  if (visualization)
178  {
179 
180  const double src_err = src_fun.ComputeL2Error(coeff);
181  const double dest_err = dest_fun.ComputeL2Error(coeff);
182 
183  if (rank == 0)
184  {
185  mfem::out << "l2 error: src: " << src_err << ", dest: " << dest_err
186  << std::endl;
187  }
188 
189  plot(*p_src_mesh, src_fun, "source");
190  plot(*p_dest_mesh, dest_fun, "destination");
191  }
192  }
193  else
194  {
195  mfem::out << "Transfer failed! Use --verbose option for diagnostic!" <<
196  std::endl;
197  }
198 
199  // Finalize transfer library context
201  return MPI_Finalize();
202 }
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
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 SetMaxSolverIterations(const int max_solver_iterations)
Control the maximum numbers of conjugate gradients steps for mass matrix inversion.
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)
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void make_fun(mfem::FiniteElementSpace &fe, mfem::Coefficient &c, mfem::GridFunction &f)
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()