libMesh
meshdiff.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 two solution files named on standard input, find all
20 // variables they have in common, and output the Hilbert norms of the
21 // differences between them.
22 
23 #include "libmesh/libmesh.h"
24 
25 #include "libmesh/mesh.h"
26 #include "libmesh/equation_systems.h"
27 #include "libmesh/exact_solution.h"
28 #include "libmesh/mesh_function.h"
29 #include "libmesh/namebased_io.h"
30 #include "libmesh/numeric_vector.h"
31 
32 using namespace libMesh;
33 
34 unsigned char dim = 2; // This gets overridden by most mesh formats
35 
36 int main(int argc, char ** argv)
37 {
38  LibMeshInit init(argc, argv);
39 
40  Mesh coarse_mesh(init.comm(), dim), fine_mesh(init.comm(), dim);
41  EquationSystems coarse_es(coarse_mesh), fine_es(fine_mesh);
42 
43  libMesh::out << "Usage: " << argv[0]
44  << " coarsemesh coarsesolution finemesh finesolution [outputdiff]" << std::endl;
45 
46  if (argc < 5)
47  libmesh_error();
48 
49  coarse_mesh.read(argv[1]);
50  libMesh::out << "Loaded coarse mesh " << argv[1] << std::endl;
51  coarse_es.read(argv[2]);
52  libMesh::out << "Loaded coarse solution " << argv[2] << std::endl;
53  fine_mesh.read(argv[3]);
54  libMesh::out << "Loaded fine mesh " << argv[3] << std::endl;
55  fine_es.read(argv[4]);
56  libMesh::out << "Loaded fine solution " << argv[4] << std::endl;
57 
58  ExactSolution exact_sol(coarse_es);
59  exact_sol.attach_reference_solution(&fine_es);
60 
61  std::vector<std::string> sysnames;
62  sysnames.reserve(coarse_es.n_systems());
63 
64  for (unsigned int i = 0; i != coarse_es.n_systems(); ++i)
65  if (fine_es.has_system(coarse_es.get_system(i).name()))
66  sysnames.push_back(coarse_es.get_system(i).name());
67  else
68  libMesh::out << "Coarse solution system "
69  << coarse_es.get_system(i).name()
70  << " not found in fine solution!" << std::endl;
71 
72  for (unsigned int i = 0; i != fine_es.n_systems(); ++i)
73  if (!coarse_es.has_system(fine_es.get_system(i).name()))
74  libMesh::out << "Fine solution system "
75  << fine_es.get_system(i).name()
76  << " not found in coarse solution!" << std::endl;
77 
78  if (!coarse_es.n_systems() && !fine_es.n_systems())
79  libMesh::out << "No systems found in fine or coarse solution!"
80  << std::endl;
81 
82  for (std::size_t i = 0; i != sysnames.size(); ++i)
83  {
84  const std::string sysname = sysnames[i];
85  const System & coarse_sys = coarse_es.get_system(sysname);
86  const System & fine_sys = fine_es.get_system(sysname);
87 
88  for (unsigned int j = 0; j != coarse_sys.n_vars(); ++j)
89  {
90  const std::string varname = coarse_sys.variable_name(j);
91 
92  if (!fine_sys.has_variable(varname))
93  {
94  libMesh::out << "Fine solution system " << sysname
95  << " variable " << varname
96  << " not found in coarse solution!" << std::endl;
97  continue;
98  }
99  const unsigned int j2 = fine_sys.variable_number(varname);
100 
101  exact_sol.compute_error(sysname, varname);
102 
103  libMesh::out << "Errors in system " << sysname << ", variable " << varname << ":" << std::endl;
104  libMesh::out << "L2 error: " << exact_sol.l2_error(sysname, varname)
105  << ", fine norm: " << coarse_sys.calculate_norm(*coarse_sys.solution, j, L2)
106  << ", coarse norm: " << fine_sys.calculate_norm(*fine_sys.solution, j2, L2) << std::endl;
107  libMesh::out << "H1 error: " << exact_sol.h1_error(sysname, varname)
108  << ", fine norm: " << coarse_sys.calculate_norm(*coarse_sys.solution, j, H1)
109  << ", coarse norm: " << fine_sys.calculate_norm(*fine_sys.solution, j2, H1) << std::endl;
110  libMesh::out << "H2 error: " << exact_sol.h2_error(sysname, varname)
111  << ", fine norm: " << coarse_sys.calculate_norm(*coarse_sys.solution, j, H2)
112  << ", coarse norm: " << fine_sys.calculate_norm(*fine_sys.solution, j2, H2) << std::endl;
113  }
114 
115  for (unsigned int j = 0; j != fine_sys.n_vars(); ++j)
116  {
117  const std::string varname = fine_sys.variable_name(j);
118 
119  if (!coarse_sys.has_variable(varname))
120  {
121  libMesh::out << "Coarse solution system " << sysname
122  << " variable " << varname
123  << " not found in fine solution!" << std::endl;
124  continue;
125  }
126  }
127 
128  if (!coarse_sys.n_vars() && !fine_sys.n_vars())
129  libMesh::out << "No variables found in fine or coarse solution system "
130  << sysname << '!' << std::endl;
131  }
132 
133  if (argc > 5)
134  {
135  libMesh::out << "Writing diff solution " << argv[5] << std::endl;
136 
137  for (std::size_t i = 0; i != sysnames.size(); ++i)
138  {
139  const std::string sysname = sysnames[i];
140  const System & coarse_sys = coarse_es.get_system(sysname);
141  const System & fine_sys = fine_es.get_system(sysname);
142 
143  UniquePtr<NumericVector<Number>> fine_solution = fine_sys.solution->clone();
144  fine_sys.solution->zero();
145 
146  std::vector<unsigned int>
147  var_remapping(fine_sys.n_vars(), libMesh::invalid_uint);
148 
149  for (unsigned int j = 0; j != fine_sys.n_vars(); ++j)
150  {
151  const std::string varname = fine_sys.variable_name(j);
152 
153  if (coarse_sys.has_variable(varname))
154  var_remapping[j] = coarse_sys.variable_number(varname);
155  }
156 
157  MeshFunction coarse_solution
158  (coarse_es, *coarse_sys.solution,
159  coarse_sys.get_dof_map(), var_remapping);
160  coarse_solution.init();
161  const DenseVector<Number> oom_value(fine_sys.n_vars(), 0);
162  coarse_solution.enable_out_of_mesh_mode(oom_value);
163 
164  fine_sys.project_solution(&coarse_solution);
165  *fine_sys.solution -= *fine_solution;
166  }
167 
168  NameBasedIO(fine_mesh).write_equation_systems (argv[5], fine_es);
169  }
170 
171  return 0;
172 }
This class handles the computation of the L2 and/or H1 error for the Systems in the EquationSystems o...
This class supports simple reads and writes in any libMesh-supported format, by dispatching to one of...
Definition: namebased_io.h:44
This is the EquationSystems class.
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2134
bool has_variable(const std::string &var) const
Definition: system.C:1256
Real h2_error(const std::string &sys_name, const std::string &unknown_name)
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr) const
Projects arbitrary functions onto the current solution.
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.
void compute_error(const std::string &sys_name, const std::string &unknown_name)
Computes and stores the error in the solution value e = u-u_h, the gradient grad(e) = grad(u) - grad(...
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
const DofMap & get_dof_map() const
Definition: system.h:2030
int main(int argc, char **argv)
Definition: meshdiff.C:36
unsigned char dim
Definition: meshdiff.C:34
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
virtual void write_equation_systems(const std::string &filename, const EquationSystems &es, const std::set< std::string > *system_names=libmesh_nullptr) libmesh_override
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: namebased_io.C:490
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
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=libmesh_nullptr) const
Definition: system.C:1383
Real h1_error(const std::string &sys_name, const std::string &unknown_name)
OStreamProxy out
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
void attach_reference_solution(const EquationSystems *es_fine)
Attach function similar to system.h which allows the user to attach a second EquationSystems object w...
Real l2_error(const std::string &sys_name, const std::string &unknown_name)
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
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