libMesh
projection.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // Open the input mesh and corresponding solution file named in command line
20 // arguments, open the output mesh, project that solution onto the
21 // output mesh, and write a corresponding output solution file.
22 
23 #include <map>
24 #include <string>
25 
26 #include "libmesh/libmesh.h"
27 
28 #include "libmesh/dof_map.h"
29 #include "libmesh/equation_systems.h"
30 #include "libmesh/getpot.h"
31 #include "libmesh/mesh.h"
32 #include "libmesh/mesh_function.h"
33 #include "libmesh/numeric_vector.h"
34 #include "libmesh/point.h"
35 #include "libmesh/replicated_mesh.h"
36 
37 
38 using namespace libMesh;
39 
40 
41 // If there's a missing input argument, then print a help message
42 void usage_error(const char * progname)
43 {
44  libMesh::out << "Options: " << progname << '\n'
45  << " --dim d mesh dimension [default: autodetect]\n"
46  << " --inmesh filename input mesh file\n"
47  << " --insoln filename input solution file\n"
48  << " --outmesh filename output mesh file [default: out_<inmesh>]\n"
49  << " --outsoln filename output solution file [default: out_<insoln>]\n"
50  << std::endl;
51 
52  exit(1);
53 }
54 
55 // Get an input argument, or print a help message if it's missing
56 template <typename T>
57 T assert_argument (GetPot & cl,
58  const std::string & argname,
59  const char * progname,
60  const T & defaultarg)
61 {
62  if (!cl.search(argname))
63  {
64  libMesh::err << ("No " + argname + " argument found!") << std::endl;
65  usage_error(progname);
66  }
67  return cl.next(defaultarg);
68 }
69 
70 // Global collections of MeshFunctions are necessary to work with
71 // global functions fptr and gptr.
72 // TODO: libMesh needs functor-based alternatives to these types of
73 // function arguments
74 std::string current_sys_name;
75 std::map<std::string, MeshFunction *> mesh_functions;
76 
77 // Return the function value on the old mesh and solution
78 Number fptr(const Point & p,
79  const Parameters &,
80  const std::string & libmesh_dbg_var(sys_name),
81  const std::string & unknown_name)
82 {
83  libmesh_assert_equal_to (sys_name, current_sys_name);
84  libmesh_assert(mesh_functions.count(unknown_name));
85  libmesh_assert(mesh_functions[unknown_name]);
86 
87  MeshFunction & meshfunc = *mesh_functions[unknown_name];
88 
89  return meshfunc(p);
90 }
91 
92 // Return the function gradient on the old mesh and solution
93 Gradient gptr(const Point & p,
94  const Parameters &,
95  const std::string & libmesh_dbg_var(sys_name),
96  const std::string & unknown_name)
97 {
98  libmesh_assert_equal_to (sys_name, current_sys_name);
99  libmesh_assert(mesh_functions.count(unknown_name));
100  libmesh_assert(mesh_functions[unknown_name]);
101 
102  MeshFunction & meshfunc = *mesh_functions[unknown_name];
103 
104  return meshfunc.gradient(p);
105 }
106 
107 
108 int main(int argc, char ** argv)
109 {
110  LibMeshInit init(argc, argv);
111 
112  GetPot cl(argc, argv);
113 
114  // In case the mesh file doesn't let us auto-infer dimension, we let
115  // the user specify it on the command line
116  const unsigned char requested_dim =
117  cast_int<unsigned char>(cl.follow(3, "--dim"));
118 
119  // Load the old mesh from --inmesh filename.
120  // Keep it serialized; we don't want elements on the new mesh to be
121  // looking for data on old mesh elements that live off-processor.
122  ReplicatedMesh old_mesh(init.comm(), requested_dim);
123 
124  const std::string meshname =
125  assert_argument(cl, "--inmesh", argv[0], std::string("mesh.xda"));
126 
127  old_mesh.read(meshname);
128  std::cout << "Old Mesh:" << std::endl;
129  old_mesh.print_info();
130 
131  // Load the new mesh from --outmesh filename
132  Mesh new_mesh(init.comm(), requested_dim);
133 
134  const std::string outmeshname = cl.follow(std::string("out_"+meshname), "--outmesh");
135 
136  new_mesh.read(outmeshname);
137  std::cout << "New Mesh:" << std::endl;
138  new_mesh.print_info();
139 
140  // Load the old solution from --insoln filename
141  // Construct the new solution from the old solution's headers, so
142  // that we get the system names and types, variable names and types,
143  // etc.
144  const std::string solnname =
145  assert_argument(cl, "--insoln", argv[0], std::string("soln.xda"));
146 
147  EquationSystems old_es(old_mesh);
148  EquationSystems new_es(new_mesh);
149 
150  XdrMODE read_mode;
151 
152  if (solnname.rfind(".xdr") < solnname.size())
153  read_mode = DECODE;
154  else if (solnname.rfind(".xda") < solnname.size())
155  read_mode = READ;
156  else
157  libmesh_error_msg("Unrecognized file extension on " << solnname);
158 
159  old_es.read(solnname, read_mode,
164 
165  new_es.read(solnname, read_mode,
168 
169  old_es.print_info();
170 
171  unsigned int n_systems = old_es.n_systems();
172  libmesh_assert_equal_to (new_es.n_systems(), n_systems);
173 
174  // For each system, serialize the solution so we can project it onto
175  // a potentially-very-different partitioning
176  for (unsigned int i = 0; i != n_systems; ++i)
177  {
178  System & old_sys = old_es.get_system(i);
179  current_sys_name = old_sys.name();
180 
181  libMesh::out << "Projecting system " << current_sys_name << std::endl;
182 
184 
185  System & new_sys = new_es.get_system(current_sys_name);
186  unsigned int n_vars = old_sys.n_vars();
187  libmesh_assert_equal_to (new_sys.n_vars(), n_vars);
188 
189  UniquePtr<NumericVector<Number>> comparison_soln =
191  std::vector<Number> global_soln;
192  old_sys.update_global_solution(global_soln);
193  comparison_soln->init(old_sys.solution->size(), true, SERIAL);
194  (*comparison_soln) = global_soln;
195 
196  // For each variable, construct a MeshFunction returning that
197  // variable's value
198  for (unsigned int j = 0; j != n_vars; ++j)
199  {
200  libMesh::out << " with variable " << old_sys.variable_name(j) << std::endl;
201 
202  MeshFunction * mesh_func =
203  new MeshFunction(old_es, *comparison_soln,
204  old_sys.get_dof_map(), j);
205  mesh_func->init();
206  mesh_functions[old_sys.variable_name(j)] = mesh_func;
207  }
208 
209  // Project all variables to the new system
210  new_sys.project_solution(fptr, gptr, old_es.parameters);
211 
212  // Clean up the MeshFunctions here so we don't bloat memory
213  for (unsigned int j = 0; j != n_vars; ++j)
214  delete mesh_functions[old_sys.variable_name(j)];
215  }
216 
217  // Write out the new solution file
218  const std::string outsolnname = cl.follow(std::string("out_"+solnname), "--outsoln");
219 
220  new_es.write(outsolnname.c_str(),
223  libMesh::out << "Wrote solution " << outsolnname << std::endl;
224 
225  return 0;
226 }
OStreamProxy err
This is the EquationSystems class.
std::map< std::string, MeshFunction * > mesh_functions
Definition: projection.C:75
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false) libmesh_override
Reads the file specified by name.
bool has_system(const std::string &name) const
void write(const std::string &name, const XdrMODE, const unsigned int write_flags=(WRITE_DATA), bool partition_agnostic=true) const
Write the systems to disk using the XDR data format.
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2134
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
Gradient gradient(const Point &p, const Real time=0.)
void read(const std::string &name, const XdrMODE, const unsigned int read_flags=(READ_HEADER|READ_DATA), bool partition_agnostic=true)
Read & initialize the systems from disk using the XDR data format.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
const std::string & name() const
Definition: system.h:1998
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
int main(int argc, char **argv)
Definition: projection.C:108
void usage_error(const char *progname)
Definition: projection.C:42
void update_global_solution(std::vector< Number > &global_soln) const
Fill the input vector global_soln so that it contains the global solution on all processors.
Definition: system.C:662
const DofMap & get_dof_map() const
Definition: system.h:2030
XdrMODE
Defines an enum for read/write mode in Xdr format.
Definition: enum_xdr_mode.h:32
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
PetscErrorCode Vec Mat libmesh_dbg_var(j)
unsigned int n_systems() const
T assert_argument(GetPot &cl, const std::string &argname, const char *progname, const T &defaultarg)
Definition: projection.C:57
std::string current_sys_name
Definition: projection.C:74
const Parallel::Communicator & comm() const
OStreamProxy out
Gradient gptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:93
Parameters parameters
Data structure holding arbitrary parameters.
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
This class provides function-like objects for data distributed over a mesh.
Definition: mesh_function.h:53
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void init() libmesh_override
Override the FunctionBase::init() member function by calling our own and specifying the Trees::NODES ...
Definition: mesh_function.h:96