libMesh
calculator.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, parse the function specified in a command line argument,
21 // L2-project its value onto the mesh, and output the new solution.
22 
23 #include <map>
24 #include <string>
25 
26 #include "L2system.h"
27 
28 #include "libmesh/libmesh.h"
29 
30 #include "libmesh/dof_map.h"
31 #include "libmesh/equation_systems.h"
32 #include "libmesh/getpot.h"
33 #include "libmesh/mesh.h"
34 #include "libmesh/newton_solver.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/parsed_fem_function.h"
37 #include "libmesh/point.h"
38 #include "libmesh/steady_solver.h"
39 
40 
41 using namespace libMesh;
42 
43 
44 // If there's a missing input argument, then print a help message
45 void usage_error(const char * progname)
46 {
47  libMesh::out << "Options: " << progname << '\n'
48  << " --dim d mesh dimension [default: autodetect]\n"
49  << " --inmesh filename input mesh file\n"
50  << " --insoln filename input solution file\n"
51  << " --calc func function to calculate\n"
52  << " --insys num input system number [default: 0]\n"
53  << " --outsoln filename output solution file [default: out_<insoln>]\n"
54  << " --family famname output FEM family [default: LAGRANGE]\n"
55  << " --order num output FEM order [default: 1]\n"
56  << std::endl;
57 
58  exit(1);
59 }
60 
61 // Get an input argument, or print a help message if it's missing
62 template <typename T>
63 T assert_argument (GetPot & cl,
64  const std::string & argname,
65  const char * progname,
66  const T & defaultarg)
67 {
68  if (!cl.search(argname))
69  {
70  libMesh::err << ("No " + argname + " argument found!") << std::endl;
71  usage_error(progname);
72  }
73  return cl.next(defaultarg);
74 }
75 
76 
77 int main(int argc, char ** argv)
78 {
79  LibMeshInit init(argc, argv);
80 
81  GetPot cl(argc, argv);
82 
83  // In case the mesh file doesn't let us auto-infer dimension, we let
84  // the user specify it on the command line
85  const unsigned char requested_dim =
86  cast_int<unsigned char>(cl.follow(3, "--dim"));
87 
88  // Load the old mesh from --inmesh filename.
89  // Keep it serialized; we don't want elements on the new mesh to be
90  // looking for data on old mesh elements that live off-processor.
91  Mesh old_mesh(init.comm(), requested_dim);
92 
93  const std::string meshname =
94  assert_argument(cl, "--inmesh", argv[0], std::string("mesh.xda"));
95 
96  old_mesh.read(meshname);
97  std::cout << "Mesh:" << std::endl;
98  old_mesh.print_info();
99 
100  // Create a new mesh for a new EquationSystems
101  Mesh new_mesh(init.comm(), requested_dim);
102  new_mesh.read(meshname);
103 
104  // Load the old solution from --insoln filename
105  // Construct the new solution from the old solution's headers, so
106  // that we get the system names and types, variable names and types,
107  // etc.
108  const std::string solnname =
109  assert_argument(cl, "--insoln", argv[0], std::string("soln.xda"));
110 
111  EquationSystems old_es(old_mesh);
112  EquationSystems new_es(new_mesh);
113 
114  XdrMODE read_mode;
115 
116  if (solnname.rfind(".xdr") < solnname.size())
117  read_mode = DECODE;
118  else if (solnname.rfind(".xda") < solnname.size())
119  read_mode = READ;
120  else
121  libmesh_error_msg("Unrecognized file extension on " << solnname);
122 
123  old_es.read(solnname, read_mode,
128 
129  old_es.print_info();
130 
131 
132  const unsigned int sysnum =
133  cl.follow(0, "--insys");
134 
135  libmesh_assert_less(sysnum, old_es.n_systems());
136 
137  System & old_sys = old_es.get_system(sysnum);
138  std::string current_sys_name = old_sys.name();
139 
140  libMesh::out << "Calculating with system " << current_sys_name << std::endl;
141 
142  L2System & new_sys = new_es.add_system<L2System>(current_sys_name);
143 
144  new_sys.time_solver =
145  UniquePtr<TimeSolver>(new SteadySolver(new_sys));
146 
147  new_sys.fe_family() =
148  cl.follow(std::string("LAGRANGE"), "--family");
149 
150  new_sys.fe_order() =
151  cl.follow(1, "--order");
152 
153  const std::string calcfunc =
154  assert_argument(cl, "--calc", argv[0], std::string(""));
155 
156  ParsedFEMFunction<Number> goal_function(old_sys, calcfunc);
157 
158  new_sys.goal_func.reset(goal_function.clone().release());
159  new_sys.input_system = &old_sys;
160 
161  new_es.init();
162 
163  DiffSolver & solver = *(new_sys.time_solver->diff_solver().get());
164  solver.quiet = false;
165  solver.verbose = true;
166  solver.relative_step_tolerance = 1e-10;
167 
168  new_sys.solve();
169 
170  // Write out the new solution file
171  const std::string outsolnname = cl.follow(std::string("out_"+solnname), "--outsoln");
172 
173  new_es.write(outsolnname.c_str(),
176  libMesh::out << "Wrote solution " << outsolnname << std::endl;
177 
178  return 0;
179 }
virtual UniquePtr< FEMFunctionBase< Output > > clone() const libmesh_override
OStreamProxy err
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
This is the EquationSystems class.
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.
unsigned int & fe_order()
Definition: L2system.h:25
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.
libMesh::System * input_system
Definition: L2system.h:34
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:162
std::string & fe_family()
Definition: L2system.h:24
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
The libMesh namespace provides an interface to certain functionality in the library.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
PetscDiffSolver & solver
libMesh::UniquePtr< libMesh::FEMFunctionBase< libMesh::Number > > goal_func
Definition: L2system.h:32
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:70
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: fem_system.C:1056
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.
ParsedFEMFunction provides support for FParser-based parsed functions in FEMSystem.
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.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
T assert_argument(GetPot &cl, const std::string &argname, const char *progname, const T &defaultarg)
Definition: calculator.C:63
unsigned int n_systems() const
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:168
std::string current_sys_name
Definition: projection.C:74
OStreamProxy out
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
int main(int argc, char **argv)
Definition: calculator.C:77
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void usage_error(const char *progname)
Definition: calculator.C:45