MFEM  v4.6.0
Finite element discretization library
field-interp.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // Field Interp Miniapp: Transfer a grid function between meshes
14 // -------------------------------------------------------------
15 //
16 // This miniapp provides the capability to transfer a grid function (H1, L2,
17 // H(div), and H(curl)) from one mesh onto another using GSLIB-FindPoints. Using
18 // FindPoints, we identify the nodal positions of the target mesh with respect
19 // to the source mesh and then interpolate the source grid function. The
20 // interpolated values are then projected onto the desired finite element space
21 // on the target mesh. Finally, the transferred solution is visualized using
22 // GLVis. Note that the source grid function can be a user-defined vector
23 // function or a grid function file that is compatible with the source mesh.
24 //
25 // Compile with: make field-interp
26 //
27 // Sample runs:
28 // field-interp
29 // field-interp -o 1
30 // field-interp -fts 3 -ft 0
31 // field-interp -m1 triple-pt-1.mesh -s1 triple-pt-1.gf -m2 triple-pt-2.mesh -ft 1
32 // field-interp -m1 triple-pt-1.mesh -m2 triple-pt-2.mesh -ft 1
33 // field-interp -m2 ../meshing/amr-quad-q2.mesh -ft 0 -r 1 -fts 0
34 
35 #include "mfem.hpp"
36 #include <fstream>
37 
38 using namespace mfem;
39 using namespace std;
40 
41 // Scalar function to project
42 double scalar_func(const Vector &x)
43 {
44  const int dim = x.Size();
45  double res = 0.0;
46  for (int d = 0; d < dim; d++) { res += x(d) * x(d); }
47  return res;
48 }
49 
50 void vector_func(const Vector &p, Vector &F)
51 {
52  F(0) = scalar_func(p);
53  for (int i = 1; i < F.Size(); i++) { F(i) = (i+1)*pow(-1, i)*F(0); }
54 }
55 
56 int main (int argc, char *argv[])
57 {
58  // Set the method's default parameters.
59  const char *src_mesh_file = "../meshing/square01.mesh";
60  const char *tar_mesh_file = "../../data/inline-tri.mesh";
61  const char *src_sltn_file = "must_be_provided_by_the_user.gf";
62  int src_fieldtype = 0;
63  int src_ncomp = 1;
64  int src_gf_ordering = 0;
65  int ref_levels = 0;
66  int fieldtype = -1;
67  int order = 3;
68  bool visualization = true;
69 
70  // Parse command-line options.
71  OptionsParser args(argc, argv);
72  args.AddOption(&src_mesh_file, "-m1", "--mesh1",
73  "Mesh file for the starting solution.");
74  args.AddOption(&tar_mesh_file, "-m2", "--mesh2",
75  "Mesh file for interpolation.");
76  args.AddOption(&src_sltn_file, "-s1", "--solution1",
77  "(optional) GridFunction file compatible with src_mesh_file."
78  "Set src_fieldtype to -1 if this option is used.");
79  args.AddOption(&src_fieldtype, "-fts", "--field-type-src",
80  "Source GridFunction type:"
81  "0 - H1 (default), 1 - L2, 2 - H(div), 3 - H(curl).");
82  args.AddOption(&src_ncomp, "-nc", "--ncomp",
83  "Number of components for H1 or L2 GridFunctions.");
84  args.AddOption(&src_gf_ordering, "-gfo", "--gfo",
85  "Node ordering: 0 (byNodes) or 1 (byVDim)");
86  args.AddOption(&ref_levels, "-r", "--refine",
87  "Number of refinements of the interpolation mesh.");
88  args.AddOption(&fieldtype, "-ft", "--field-type",
89  "Target GridFunction type: -1 - source GridFunction type (default),"
90  "0 - H1, 1 - L2, 2 - H(div), 3 - H(curl).");
91  args.AddOption(&order, "-o", "--order",
92  "Order of the interpolated solution.");
93  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
94  "--no-visualization",
95  "Enable or disable GLVis visualization.");
96  args.Parse();
97  if (!args.Good())
98  {
99  args.PrintUsage(cout);
100  return 1;
101  }
102  args.PrintOptions(cout);
103 
104  // If a gridfunction is specified, set src_fieldtype to -1
105  if (strcmp(src_sltn_file, "must_be_provided_by_the_user.gf") != 0)
106  {
107  src_fieldtype = -1;
108  }
109 
110  // Input meshes.
111  Mesh mesh_1(src_mesh_file, 1, 1, false);
112  Mesh mesh_2(tar_mesh_file, 1, 1, false);
113  const int dim = mesh_1.Dimension();
114  MFEM_ASSERT(dim == mesh_2.Dimension(), "Source and target meshes "
115  "must be in the same dimension.");
116  MFEM_VERIFY(dim > 1, "GSLIB requires a 2D or a 3D mesh" );
117 
118  for (int lev = 0; lev < ref_levels; lev++)
119  {
120  mesh_2.UniformRefinement();
121  }
122 
123  if (mesh_1.GetNodes() == NULL) { mesh_1.SetCurvature(1); }
124  if (mesh_2.GetNodes() == NULL) { mesh_2.SetCurvature(1); }
125  const int mesh_poly_deg =
126  mesh_2.GetNodes()->FESpace()->GetElementOrder(0);
127  cout << "Source mesh curvature: "
128  << mesh_1.GetNodes()->OwnFEC()->Name() << endl
129  << "Target mesh curvature: "
130  << mesh_2.GetNodes()->OwnFEC()->Name() << endl;
131 
132  int src_vdim = src_ncomp;
133  FiniteElementCollection *src_fec = NULL;
134  FiniteElementSpace *src_fes = NULL;
135  GridFunction *func_source = NULL;
136  if (src_fieldtype < 0) // use src_sltn_file
137  {
138  ifstream mat_stream_1(src_sltn_file);
139  func_source = new GridFunction(&mesh_1, mat_stream_1);
140  src_vdim = func_source->FESpace()->GetVDim();
141  src_fes = func_source->FESpace();
142  }
143  else if (src_fieldtype == 0)
144  {
145  src_fec = new H1_FECollection(order, dim);
146  }
147  else if (src_fieldtype == 1)
148  {
149  src_fec = new L2_FECollection(order, dim);
150  }
151  else if (src_fieldtype == 2)
152  {
153  src_fec = new RT_FECollection(order, dim);
154  src_ncomp = 1;
155  src_vdim = dim;
156  }
157  else if (src_fieldtype == 3)
158  {
159  src_fec = new ND_FECollection(order, dim);
160  src_ncomp = 1;
161  src_vdim = dim;
162  }
163  else
164  {
165  MFEM_ABORT("Invalid FECollection type.");
166  }
167 
168  if (src_fieldtype > -1)
169  {
170  MFEM_VERIFY(src_gf_ordering >= 0 && src_gf_ordering <= 1,
171  " Source grid function ordering must be 0 (byNodes)"
172  " or 1 (byVDIM.");
173  src_fes = new FiniteElementSpace(&mesh_1, src_fec, src_ncomp, src_gf_ordering);
174  func_source = new GridFunction(src_fes);
175  // Project the grid function using VectorFunctionCoefficient.
177  func_source->ProjectCoefficient(F);
178  }
179 
180  // Display the starting mesh and the field.
181  if (visualization)
182  {
183  char vishost[] = "localhost";
184  int visport = 19916;
185  socketstream sout1;
186  sout1.open(vishost, visport);
187  if (!sout1)
188  {
189  cout << "Unable to connect to GLVis server at "
190  << vishost << ':' << visport << endl;
191  }
192  else
193  {
194  sout1.precision(8);
195  sout1 << "solution\n" << mesh_1 << *func_source
196  << "window_title 'Source mesh and solution'"
197  << "window_geometry 0 0 600 600";
198  if (dim == 2) { sout1 << "keys RmjAc"; }
199  if (dim == 3) { sout1 << "keys mA\n"; }
200  sout1 << flush;
201  }
202  }
203 
204  const Geometry::Type gt = mesh_2.GetNodalFESpace()->GetFE(0)->GetGeomType();
205  MFEM_VERIFY(gt != Geometry::PRISM, "Wedge elements are not currently "
206  "supported.");
207  MFEM_VERIFY(mesh_2.GetNumGeometries(mesh_2.Dimension()) == 1, "Mixed meshes"
208  "are not currently supported.");
209 
210  // Ensure the source grid function can be transferred using GSLIB-FindPoints.
211  const FiniteElementCollection *fec_in = func_source->FESpace()->FEColl();
212  std::cout << "Source FE collection: " << fec_in->Name() << std::endl;
213 
214  if (src_fieldtype < 0)
215  {
216  const H1_FECollection *fec_h1 = dynamic_cast<const H1_FECollection *>(fec_in);
217  const L2_FECollection *fec_l2 = dynamic_cast<const L2_FECollection *>(fec_in);
218  const RT_FECollection *fec_rt = dynamic_cast<const RT_FECollection *>(fec_in);
219  const ND_FECollection *fec_nd = dynamic_cast<const ND_FECollection *>(fec_in);
220  if (fec_h1) { src_fieldtype = 0; }
221  else if (fec_l2) { src_fieldtype = 1; }
222  else if (fec_rt) { src_fieldtype = 2; }
223  else if (fec_nd) { src_fieldtype = 3; }
224  else { MFEM_ABORT("GridFunction type not supported yet."); }
225  }
226  if (fieldtype < 0) { fieldtype = src_fieldtype; }
227 
228  // Setup the FiniteElementSpace and GridFunction on the target mesh.
229  FiniteElementCollection *tar_fec = NULL;
230  FiniteElementSpace *tar_fes = NULL;
231 
232  int tar_vdim = src_vdim;
233  if (fieldtype == 0)
234  {
235  tar_fec = new H1_FECollection(order, dim);
236  tar_vdim = (src_fieldtype > 1) ? dim : src_vdim;
237  }
238  else if (fieldtype == 1)
239  {
240  tar_fec = new L2_FECollection(order, dim);
241  tar_vdim = (src_fieldtype > 1) ? dim : src_vdim;
242  }
243  else if (fieldtype == 2)
244  {
245  tar_fec = new RT_FECollection(order, dim);
246  tar_vdim = 1;
247  MFEM_VERIFY(src_fieldtype > 1, "Cannot interpolate a scalar "
248  "grid function to a vector");
249 
250  }
251  else if (fieldtype == 3)
252  {
253  tar_fec = new ND_FECollection(order, dim);
254  tar_vdim = 1;
255  MFEM_VERIFY(src_fieldtype > 1, "Cannot interpolate a scalar "
256  "grid function to a vector");
257  }
258  else
259  {
260  MFEM_ABORT("GridFunction type not supported.");
261  }
262  std::cout << "Target FE collection: " << tar_fec->Name() << std::endl;
263  tar_fes = new FiniteElementSpace(&mesh_2, tar_fec, tar_vdim,
264  src_fes->GetOrdering());
265  GridFunction func_target(tar_fes);
266 
267  const int NE = mesh_2.GetNE(),
268  nsp = tar_fes->GetFE(0)->GetNodes().GetNPoints(),
269  tar_ncomp = func_target.VectorDim();
270 
271  // Generate list of points where the grid function will be evaluated.
272  Vector vxyz;
273  int point_ordering;
274  if (fieldtype == 0 && order == mesh_poly_deg)
275  {
276  vxyz = *mesh_2.GetNodes();
277  point_ordering = mesh_2.GetNodes()->FESpace()->GetOrdering();
278  }
279  else
280  {
281  vxyz.SetSize(nsp*NE*dim);
282  for (int i = 0; i < NE; i++)
283  {
284  const FiniteElement *fe = tar_fes->GetFE(i);
285  const IntegrationRule ir = fe->GetNodes();
287 
288  DenseMatrix pos;
289  et->Transform(ir, pos);
290  Vector rowx(vxyz.GetData() + i*nsp, nsp),
291  rowy(vxyz.GetData() + i*nsp + NE*nsp, nsp),
292  rowz;
293  if (dim == 3)
294  {
295  rowz.SetDataAndSize(vxyz.GetData() + i*nsp + 2*NE*nsp, nsp);
296  }
297  pos.GetRow(0, rowx);
298  pos.GetRow(1, rowy);
299  if (dim == 3) { pos.GetRow(2, rowz); }
300  }
301  point_ordering = Ordering::byNODES;
302  }
303  const int nodes_cnt = vxyz.Size() / dim;
304 
305  // Evaluate source grid function.
306  Vector interp_vals(nodes_cnt*tar_ncomp);
307  FindPointsGSLIB finder;
308  finder.Setup(mesh_1);
309  finder.Interpolate(vxyz, *func_source, interp_vals, point_ordering);
310 
311  // Project the interpolated values to the target FiniteElementSpace.
312  if (fieldtype <= 1) // H1 or L2
313  {
314  if ((fieldtype == 0 && order == mesh_poly_deg) || fieldtype == 1)
315  {
316  func_target = interp_vals;
317  }
318  else // H1 - but mesh order != GridFunction order
319  {
320  Array<int> vdofs;
321  Vector vals;
322  Vector elem_dof_vals(nsp*tar_ncomp);
323 
324  for (int i = 0; i < mesh_2.GetNE(); i++)
325  {
326  tar_fes->GetElementVDofs(i, vdofs);
327  vals.SetSize(vdofs.Size());
328  for (int j = 0; j < nsp; j++)
329  {
330  for (int d = 0; d < tar_ncomp; d++)
331  {
332  // Arrange values byNodes
333  int idx = src_fes->GetOrdering() == Ordering::byNODES ?
334  d*nsp*NE + i*nsp + j : i*nsp*dim + d + j*dim;
335  elem_dof_vals(j + d*nsp) = interp_vals(idx);
336  }
337  }
338  func_target.SetSubVector(vdofs, elem_dof_vals);
339  }
340  }
341  }
342  else // H(div) or H(curl)
343  {
344  Array<int> vdofs;
345  Vector vals;
346  Vector elem_dof_vals(nsp*tar_ncomp);
347 
348  for (int i = 0; i < mesh_2.GetNE(); i++)
349  {
350  tar_fes->GetElementVDofs(i, vdofs);
351  vals.SetSize(vdofs.Size());
352  for (int j = 0; j < nsp; j++)
353  {
354  for (int d = 0; d < tar_ncomp; d++)
355  {
356  // Arrange values byVDim
357  int idx = src_fes->GetOrdering() == Ordering::byNODES ?
358  d*nsp*NE + i*nsp + j : i*nsp*dim + d + j*dim;
359  elem_dof_vals(j*tar_ncomp+d) = interp_vals(idx);
360  }
361  }
362  tar_fes->GetFE(i)->ProjectFromNodes(elem_dof_vals,
363  *tar_fes->GetElementTransformation(i),
364  vals);
365  func_target.SetSubVector(vdofs, vals);
366  }
367  }
368 
369  // Visualize the transferred solution.
370  if (visualization)
371  {
372  char vishost[] = "localhost";
373  int visport = 19916;
374  socketstream sout1;
375  sout1.open(vishost, visport);
376  if (!sout1)
377  {
378  cout << "Unable to connect to GLVis server at "
379  << vishost << ':' << visport << endl;
380  }
381  else
382  {
383  sout1.precision(8);
384  sout1 << "solution\n" << mesh_2 << func_target
385  << "window_title 'Target mesh and solution'"
386  << "window_geometry 600 0 600 600";
387  if (dim == 2) { sout1 << "keys RmjAc"; }
388  if (dim == 3) { sout1 << "keys mA\n"; }
389  sout1 << flush;
390  }
391  }
392 
393  // Output the target mesh with the interpolated solution.
394  ostringstream rho_name;
395  rho_name << "interpolated.gf";
396  ofstream rho_ofs(rho_name.str().c_str());
397  rho_ofs.precision(8);
398  func_target.Save(rho_ofs);
399  rho_ofs.close();
400 
401  // Free the internal gslib data.
402  finder.FreeData();
403 
404  // Delete remaining memory.
405  if (func_source->OwnFEC())
406  {
407  delete func_source;
408  }
409  else
410  {
411  delete func_source;
412  delete src_fes;
413  delete src_fec;
414  }
415  delete tar_fes;
416  delete tar_fec;
417 
418  return 0;
419 }
Abstract class for all finite elements.
Definition: fe_base.hpp:233
const char vishost[]
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:820
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:5630
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector of values at the finite element nodes and a transformation, compute its projection (ap...
Definition: fe_base.cpp:138
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:6274
int close()
Close the socketstream.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
double scalar_func(const Vector &x)
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
virtual void FreeData()
Definition: gslib.cpp:272
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
void Setup(Mesh &m, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
Definition: gslib.cpp:107
virtual const char * Name() const
Definition: fe_coll.hpp:80
int main(int argc, char *argv[])
const int visport
FiniteElementCollection * OwnFEC()
Definition: gridfunc.hpp:124
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
A general vector function coefficient.
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1314
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
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:154
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
void vector_func(const Vector &p, Vector &F)
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
int dim
Definition: ex24.cpp:53
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2415
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points...
Definition: gslib.hpp:52
Vector data type.
Definition: vector.hpp:58
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:389
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769