MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
lor-transfer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 
48 using namespace std;
49 using namespace mfem;
50 
51 int problem = 1; // problem type
52 
53 int Wx = 0, Wy = 0; // window position
54 int Ww = 350, Wh = 350; // window size
55 int offx = Ww+5, offy = Wh+25; // window offsets
56 
57 string space;
58 string direction;
59 
60 // Exact functions to project
61 double RHO_exact(const Vector &x);
62 
63 // Helper functions
64 void visualize(VisItDataCollection &, string, int, int);
66  string);
67 
68 int 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.Parse();
99  if (!args.Good())
100  {
101  args.PrintUsage(cout);
102  return 1;
103  }
104  args.PrintOptions(cout);
105 
106  // Read the mesh from the given mesh file.
107  Mesh mesh(mesh_file, 1, 1);
108  int dim = mesh.Dimension();
109 
110  // Create the low-order refined mesh
111  int basis_lor = BasisType::GaussLobatto; // BasisType::ClosedUniform;
112  Mesh mesh_lor = Mesh::MakeRefined(mesh, lref, basis_lor);
113 
114  // Create spaces
115  FiniteElementCollection *fec, *fec_lor;
116  if (useH1)
117  {
118  space = "H1";
119  if (lorder == 0)
120  {
121  lorder = 1;
122  cerr << "Switching the H1 LOR space order from 0 to 1\n";
123  }
124  fec = new H1_FECollection(order, dim);
125  fec_lor = new H1_FECollection(lorder, dim);
126  }
127  else
128  {
129  space = "L2";
130  fec = new L2_FECollection(order, dim);
131  fec_lor = new L2_FECollection(lorder, dim);
132  }
133 
134  FiniteElementSpace fespace(&mesh, fec);
135  FiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
136 
137  GridFunction rho(&fespace);
138  GridFunction rho_lor(&fespace_lor);
139 
140  // Data collections for vis/analysis
141  VisItDataCollection HO_dc("HO", &mesh);
142  HO_dc.RegisterField("density", &rho);
143  VisItDataCollection LOR_dc("LOR", &mesh_lor);
144  LOR_dc.RegisterField("density", &rho_lor);
145 
146  BilinearForm M_ho(&fespace);
148  M_ho.Assemble();
149  M_ho.Finalize();
150 
151  BilinearForm M_lor(&fespace_lor);
153  M_lor.Assemble();
154  M_lor.Finalize();
155 
156  // HO projections
157  direction = "HO -> LOR @ HO";
159  rho.ProjectCoefficient(RHO);
160  double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO ");
161  if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; }
162 
163  GridTransfer *gt;
164  if (use_pointwise_transfer)
165  {
166  gt = new InterpolationGridTransfer(fespace, fespace_lor);
167  }
168  else
169  {
170  gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
171  }
172  const Operator &R = gt->ForwardOperator();
173  const Operator &P = gt->BackwardOperator();
174 
175  // HO->LOR restriction
176  direction = "HO -> LOR @ LOR";
177  R.Mult(rho, rho_lor);
178  compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) ");
179  if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; }
180 
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  // HO* to LOR* dual fields
193  GridFunction ones(&fespace), ones_lor(&fespace_lor);
194  ones = 1.0;
195  ones_lor = 1.0;
196  LinearForm M_rho(&fespace), M_rho_lor(&fespace_lor);
197  if (!use_pointwise_transfer)
198  {
199  M_ho.Mult(rho, M_rho);
200  P.MultTranspose(M_rho, M_rho_lor);
201  cout << "HO -> LOR dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor))
202  << endl << endl;
203  }
204 
205  // LOR projections
206  direction = "LOR -> HO @ LOR";
207  rho_lor.ProjectCoefficient(RHO);
208  GridFunction rho_lor_prev = rho_lor;
209  double lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR ");
210  if (vis) { visualize(LOR_dc, "LOR", Wx, Wy); Wx += offx; }
211 
212  // Prolongate to HO space
213  direction = "LOR -> HO @ HO";
214  P.Mult(rho_lor, rho);
215  compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) ");
216  if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy); Wx += offx; }
217 
218  // Restrict back to LOR space. This won't give the original function because
219  // the rho_lor doesn't necessarily live in the range of R.
220  direction = "LOR -> HO @ LOR";
221  R.Mult(rho, rho_lor);
222  compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
223  if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); }
224 
225  rho_lor_prev -= rho_lor;
226  cout.precision(12);
227  cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl;
228 
229  // LOR* to HO* dual fields
230  if (!use_pointwise_transfer)
231  {
232  M_lor.Mult(rho_lor, M_rho_lor);
233  R.MultTranspose(M_rho_lor, M_rho);
234  cout << "LOR -> HO dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor))
235  << '\n';
236  }
237 
238  delete fec;
239  delete fec_lor;
240 
241  return 0;
242 }
243 
244 
245 double RHO_exact(const Vector &x)
246 {
247  switch (problem)
248  {
249  case 1: // smooth field
250  return x(1)+0.25*cos(2*M_PI*x.Norml2());
251  case 2: // cubic function
252  return x(1)*x(1)*x(1) + 2*x(0)*x(1) + x(0);
253  case 3: // sharp gradient
254  return M_PI/2-atan(5*(2*x.Norml2()-1));
255  case 4: // basis function
256  return (x.Norml2() < 0.1) ? 1 : 0;
257  default:
258  return 1.0;
259  }
260 }
261 
262 
263 void visualize(VisItDataCollection &dc, string prefix, int x, int y)
264 {
265  int w = Ww, h = Wh;
266 
267  char vishost[] = "localhost";
268  int visport = 19916;
269 
270  socketstream sol_sockL2(vishost, visport);
271  sol_sockL2.precision(8);
272  sol_sockL2 << "solution\n" << *dc.GetMesh() << *dc.GetField("density")
273  << "window_geometry " << x << " " << y << " " << w << " " << h
274  << "plot_caption '" << space << " " << prefix << " Density'"
275  << "window_title '" << direction << "'" << flush;
276 }
277 
278 
279 double compute_mass(FiniteElementSpace *L2, double massL2,
280  VisItDataCollection &dc, string prefix)
281 {
282  ConstantCoefficient one(1.0);
283  BilinearForm ML2(L2);
284  ML2.AddDomainIntegrator(new MassIntegrator(one));
285  ML2.Assemble();
286 
287  GridFunction rhoone(L2);
288  rhoone = 1.0;
289 
290  double newmass = ML2.InnerProduct(*dc.GetField("density"),rhoone);
291  cout.precision(18);
292  cout << space << " " << prefix << " mass = " << newmass;
293  if (massL2 >= 0)
294  {
295  cout.precision(4);
296  cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)";
297  }
298  cout << endl;
299  return newmass;
300 }
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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:78
GridFunction * GetField(const std::string &field_name)
Get a pointer to a grid function in the collection.
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
int offx
int Wx
virtual const Operator & BackwardOperator()=0
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:810
double RHO_exact(const Vector &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 Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
virtual const Operator & ForwardOperator()=0
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
constexpr int visport
Data collection with VisIt I/O routines.
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:379
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
Base class for transfer algorithms that construct transfer Operators between two finite element (FE) ...
Definition: transfer.hpp:29
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
Transfer data between a coarse mesh and an embedded refined mesh using interpolation.
Definition: transfer.hpp:121
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
int problem
Definition: ex15.cpp:62
string space
double InnerProduct(const Vector &x, const Vector &y) const
Compute .
string direction
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int Wh
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Transfer data in L2 and H1 finite element spaces between a coarse mesh and an embedded refined mesh u...
Definition: transfer.hpp:168
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2395
Mesh * GetMesh()
Get a pointer to the mesh in the collection.
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
int offy
int Wy
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
int Ww
Abstract operator.
Definition: operator.hpp:24
double compute_mass(FiniteElementSpace *, double, VisItDataCollection &, string)
FDualNumber< tbase > atan(const FDualNumber< tbase > &f)
atan([dual number])
Definition: fdual.hpp:455
int main()
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150