MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
lor-transfer.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// LOR Transfer Miniapp: Map functions between HO and LOR spaces
14// --------------------------------------------------------------
15//
16// This miniapp visualizes the maps between a high-order (HO) finite element
17// space, typically using high-order functions on a high-order mesh, and a
18// low-order refined (LOR) finite element space, typically defined by 0th or 1st
19// order functions on a low-order refinement of the HO mesh.
20//
21// The grid transfer operators are represented using either
22// InterpolationGridTransfer or L2ProjectionGridTransfer (depending on the
23// options requested by the user). The two transfer operators are then:
24//
25// 1. R: HO -> LOR, defined by GridTransfer::ForwardOperator
26// 2. P: LOR -> HO, defined by GridTransfer::BackwardOperator
27//
28// While defined generally, these operators have some nice properties for
29// particular finite element spaces. For example they satisfy PR=I, plus mass
30// conservation in both directions for L2 fields.
31//
32// Compile with: make lor-transfer
33//
34// Sample runs: lor-transfer
35// lor-transfer -h1
36// lor-transfer -t
37// lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
38// lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
39// lor-transfer -lref 4 -o 4 -lo 0 -p 1
40// lor-transfer -lref 5 -o 4 -lo 0 -p 1
41// lor-transfer -lref 5 -o 4 -lo 3 -p 2
42// lor-transfer -lref 5 -o 4 -lo 0 -p 3
43
44#include "mfem.hpp"
45#include <fstream>
46#include <iostream>
47
48using namespace std;
49using namespace mfem;
50
51int problem = 1; // problem type
52
53int Wx = 0, Wy = 0; // window position
54int Ww = 350, Wh = 350; // window size
55int offx = Ww+5, offy = Wh+25; // window offsets
56
57string space;
58string direction;
59
60// Exact functions to project
61real_t RHO_exact(const Vector &x);
62
63// Helper functions
64void visualize(VisItDataCollection &, string, int, int);
66 string);
67
68int main(int argc, char *argv[])
69{
70 // Parse command-line options.
71 const char *mesh_file = "../../data/star.mesh";
72 int order = 3;
73 int lref = order+1;
74 int lorder = 0;
75 bool vis = true;
76 bool useH1 = false;
77 bool use_pointwise_transfer = false;
78
79 OptionsParser args(argc, argv);
80 args.AddOption(&mesh_file, "-m", "--mesh",
81 "Mesh file to use.");
82 args.AddOption(&problem, "-p", "--problem",
83 "Problem type (see the RHO_exact function).");
84 args.AddOption(&order, "-o", "--order",
85 "Finite element order (polynomial degree) or -1 for"
86 " isoparametric space.");
87 args.AddOption(&lref, "-lref", "--lor-ref-level", "LOR refinement level.");
88 args.AddOption(&lorder, "-lo", "--lor-order",
89 "LOR space order (polynomial degree, zero by default).");
90 args.AddOption(&vis, "-vis", "--visualization", "-no-vis",
91 "--no-visualization",
92 "Enable or disable GLVis visualization.");
93 args.AddOption(&useH1, "-h1", "--use-h1", "-l2", "--use-l2",
94 "Use H1 spaces instead of L2.");
95 args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
96 "-no-t", "--dont-use-pointwise-transfer",
97 "Use pointwise transfer operators instead of L2 projection.");
98 args.ParseCheck();
99
100 // Read the mesh from the given mesh file.
101 Mesh mesh(mesh_file, 1, 1);
102 int dim = mesh.Dimension();
103
104 // Create the low-order refined mesh
105 int basis_lor = BasisType::GaussLobatto; // BasisType::ClosedUniform;
106 Mesh mesh_lor = Mesh::MakeRefined(mesh, lref, basis_lor);
107
108 // Create spaces
109 FiniteElementCollection *fec, *fec_lor;
110 if (useH1)
111 {
112 space = "H1";
113 if (lorder == 0)
114 {
115 lorder = 1;
116 cerr << "Switching the H1 LOR space order from 0 to 1\n";
117 }
118 fec = new H1_FECollection(order, dim);
119 fec_lor = new H1_FECollection(lorder, dim);
120 }
121 else
122 {
123 space = "L2";
124 fec = new L2_FECollection(order, dim);
125 fec_lor = new L2_FECollection(lorder, dim);
126 }
127
128 FiniteElementSpace fespace(&mesh, fec);
129 FiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
130
131 GridFunction rho(&fespace);
132 GridFunction rho_lor(&fespace_lor);
133
134 // Data collections for vis/analysis
135 VisItDataCollection HO_dc("HO", &mesh);
136 HO_dc.RegisterField("density", &rho);
137 VisItDataCollection LOR_dc("LOR", &mesh_lor);
138 LOR_dc.RegisterField("density", &rho_lor);
139
140 BilinearForm M_ho(&fespace);
142 M_ho.Assemble();
143 M_ho.Finalize();
144
145 BilinearForm M_lor(&fespace_lor);
147 M_lor.Assemble();
148 M_lor.Finalize();
149
150 // HO projections
151 direction = "HO -> LOR @ HO";
153 rho.ProjectCoefficient(RHO);
154 // Make sure AMR constraints are satisfied
155 rho.SetTrueVector();
156 rho.SetFromTrueVector();
157
158 real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO ");
159 if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; }
160
161 GridTransfer *gt;
162 if (use_pointwise_transfer)
163 {
164 gt = new InterpolationGridTransfer(fespace, fespace_lor);
165 }
166 else
167 {
168 gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
169 }
170 const Operator &R = gt->ForwardOperator();
171
172 // HO->LOR restriction
173 direction = "HO -> LOR @ LOR";
174 R.Mult(rho, rho_lor);
175 compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) ");
176 if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; }
177
179 {
180 const Operator &P = gt->BackwardOperator();
181 // LOR->HO prolongation
182 direction = "HO -> LOR @ HO";
183 GridFunction rho_prev = rho;
184 P.Mult(rho_lor, rho);
185 compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
186 if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy); Wx = 0; Wy += offy; }
187
188 rho_prev -= rho;
189 cout.precision(12);
190 cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl;
191 }
192
193 // HO* to LOR* dual fields
194 LinearForm M_rho(&fespace), M_rho_lor(&fespace_lor);
195 if (!use_pointwise_transfer && gt->SupportsBackwardsOperator())
196 {
197 const Operator &P = gt->BackwardOperator();
198 M_ho.Mult(rho, M_rho);
199 P.MultTranspose(M_rho, M_rho_lor);
200 cout << "HO -> LOR dual field: " << abs(M_rho.Sum()-M_rho_lor.Sum()) << "\n\n";
201 }
202
203 // LOR projections
204 direction = "LOR -> HO @ LOR";
205 rho_lor.ProjectCoefficient(RHO);
206 GridFunction rho_lor_prev = rho_lor;
207 real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR ");
208 if (vis) { visualize(LOR_dc, "LOR", Wx, Wy); Wx += offx; }
209
211 {
212 const Operator &P = gt->BackwardOperator();
213 // Prolongate to HO space
214 direction = "LOR -> HO @ HO";
215 P.Mult(rho_lor, rho);
216 compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) ");
217 if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy); Wx += offx; }
218
219 // Restrict back to LOR space. This won't give the original function because
220 // the rho_lor doesn't necessarily live in the range of R.
221 direction = "LOR -> HO @ LOR";
222 R.Mult(rho, rho_lor);
223 compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
224 if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); }
225
226 rho_lor_prev -= rho_lor;
227 cout.precision(12);
228 cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl;
229 }
230
231 // LOR* to HO* dual fields
232 if (!use_pointwise_transfer)
233 {
234 M_lor.Mult(rho_lor, M_rho_lor);
235 R.MultTranspose(M_rho_lor, M_rho);
236 cout << "LOR -> HO dual field: " << abs(M_rho.Sum() - M_rho_lor.Sum()) << '\n';
237 }
238
239 delete fec;
240 delete fec_lor;
241
242 return 0;
243}
244
245
247{
248 switch (problem)
249 {
250 case 1: // smooth field
251 return x(1)+0.25*cos(2*M_PI*x.Norml2());
252 case 2: // cubic function
253 return x(1)*x(1)*x(1) + 2*x(0)*x(1) + x(0);
254 case 3: // sharp gradient
255 return M_PI/2-atan(5*(2*x.Norml2()-1));
256 case 4: // basis function
257 return (x.Norml2() < 0.1) ? 1 : 0;
258 default:
259 return 1.0;
260 }
261}
262
263
264void visualize(VisItDataCollection &dc, string prefix, int x, int y)
265{
266 int w = Ww, h = Wh;
267
268 char vishost[] = "localhost";
269 int visport = 19916;
270
271 socketstream sol_sockL2(vishost, visport);
272 sol_sockL2.precision(8);
273 sol_sockL2 << "solution\n" << *dc.GetMesh() << *dc.GetField("density")
274 << "window_geometry " << x << " " << y << " " << w << " " << h
275 << "plot_caption '" << space << " " << prefix << " Density'"
276 << "window_title '" << direction << "'" << flush;
277}
278
279
281 VisItDataCollection &dc, string prefix)
282{
283 ConstantCoefficient one(1.0);
284 LinearForm lf(L2);
286 lf.Assemble();
287
288 real_t newmass = lf(*dc.GetField("density"));
289 cout.precision(18);
290 cout << space << " " << prefix << " mass = " << newmass;
291 if (massL2 >= 0)
292 {
293 cout.precision(4);
294 cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)";
295 }
296 cout << endl;
297 return newmass;
298}
@ GaussLobatto
Closed type.
Definition: fe_base.hpp:32
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:85
GridFunction * GetField(const std::string &field_name)
Get a pointer to a grid function in the collection.
Mesh * GetMesh()
Get a pointer to the mesh in the collection.
Class for domain integration .
Definition: lininteg.hpp:109
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
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
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2346
Base class for transfer algorithms that construct transfer Operators between two finite element (FE) ...
Definition: transfer.hpp:30
virtual bool SupportsBackwardsOperator() const
Definition: transfer.hpp:102
virtual const Operator & ForwardOperator()=0
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
virtual const Operator & BackwardOperator()=0
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:260
Transfer data between a coarse mesh and an embedded refined mesh using interpolation.
Definition: transfer.hpp:124
Transfer data in L2 and H1 finite element spaces between a coarse mesh and an embedded refined mesh u...
Definition: transfer.hpp:171
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:25
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:45
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:172
Mesh data type.
Definition: mesh.hpp:56
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
static Mesh MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
Create a refined (by any factor) version of orig_mesh.
Definition: mesh.cpp:4290
Abstract operator.
Definition: operator.hpp:25
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:255
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
Vector data type.
Definition: vector.hpp:80
real_t Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:850
real_t Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:832
real_t Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:1283
Data collection with VisIt I/O routines.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int dim
Definition: ex24.cpp:53
int visport
char vishost[]
int main()
int Ww
int Wy
void visualize(VisItDataCollection &, string, int, int)
int Wx
int problem
int offx
int Wh
string direction
real_t RHO_exact(const Vector &x)
string space
int offy
real_t compute_mass(FiniteElementSpace *, real_t, VisItDataCollection &, string)
float real_t
Definition: config.hpp:43