MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
23using namespace mfem;
24using namespace std;
25
26void destination_transform(const Vector &x, Vector &x_new)
27{
28 x_new = x;
29 // x_new *= 0.5;
30}
31
32int 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
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}
A general function coefficient.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
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
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void Update() override
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition pgridfunc.cpp:90
This class implements the parallel variational transfer between finite element spaces....
void SetAssembleMassAndCouplingTogether(const bool value)
Control if the Mass matrix is computed together with the coupling operator every time.
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 SetVerbose(const bool verbose)
Expose process details with verbose output.
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 vector function coefficient.
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
void plot(mfem::Mesh &mesh, mfem::GridFunction &x, std::string title)
void vector_fun(const mfem::Vector &x, mfem::Vector &f)
double example_fun(const mfem::Vector &x)
void check_options(mfem::OptionsParser &args)
void destination_transform(const Vector &x, Vector &x_new)
Definition ex1p.cpp:26
int main()
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
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 InitTransfer(int argc, char *argv[])
Initializes the par_moonolith library. It also calls MPI_Init.
Definition transfer.cpp:19
int FinalizeTransfer()
Finalize the par_moonolith library.
Definition transfer.cpp:23