libMesh
solution_components.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 mesh and solution file given, create a new solution file,
20 // and copy all listed variables from the old solution to the new.
21 
22 #include "libmesh/libmesh.h"
23 
24 #include "libmesh/equation_systems.h"
25 #include "libmesh/mesh.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/id_types.h"
28 #include "libmesh/elem.h"
29 
30 unsigned char dim = 2; // This gets overridden by most mesh formats
31 
32 int main(int argc, char ** argv)
33 {
34  using namespace libMesh;
35 
36  LibMeshInit init(argc, argv);
37 
38  Mesh mesh1(init.comm(), dim);
39  EquationSystems es1(mesh1);
40 
41  libMesh::out << "Usage: " << argv[0]
42  << " mesh oldsolution newsolution system1 variable1 [sys2 var2...]" << std::endl;
43 
44  // We should have one system name for each variable name, and those
45  // get preceded by an even number of arguments.
46  libmesh_assert (!(argc % 2));
47 
48  // We should have at least one system/variable pair following the
49  // initial arguments
50  libmesh_assert_greater_equal (argc, 6);
51 
52  mesh1.read(argv[1]);
53  libMesh::out << "Loaded mesh " << argv[1] << std::endl;
54  Mesh mesh2(mesh1);
55  EquationSystems es2(mesh2);
56 
57  es1.read(argv[2],
62  libMesh::out << "Loaded solution " << argv[2] << std::endl;
63 
64  std::vector<unsigned int> old_sys_num((argc-4)/2),
65  new_sys_num((argc-4)/2),
66  old_var_num((argc-4)/2),
67  new_var_num((argc-4)/2);
68 
69  std::vector<const System *> old_system((argc-4)/2);
70  std::vector<System *> new_system((argc-4)/2);
71 
72  for (int argi = 4; argi < argc; argi += 2)
73  {
74  const char * sysname = argv[argi];
75  const char * varname = argv[argi+1];
76 
77  const unsigned int pairnum = (argi-4)/2;
78 
79  libmesh_assert(es1.has_system(sysname));
80 
81  const System & old_sys = es1.get_system(sysname);
82  old_system[pairnum] = &old_sys;
83  old_sys_num[pairnum] = old_sys.number();
84 
85  libmesh_assert(old_sys.has_variable(varname));
86 
87  old_var_num[pairnum] = old_sys.variable_number(varname);
88 
89  const Variable & variable = old_sys.variable(old_var_num[pairnum]);
90 
91  std::string systype = old_sys.system_type();
92 
93  System & new_sys = es2.add_system(systype, sysname);
94  new_system[pairnum] = &new_sys;
95  new_sys_num[pairnum] = new_sys.number();
96 
97  new_var_num[pairnum] =
98  new_sys.add_variable(varname, variable.type(),
99  &variable.active_subdomains());
100  }
101 
102  es2.init();
103 
104  // A future version of this app should copy variables for
105  // non-solution vectors too.
106 
107  // Copy over any nodal degree of freedom coefficients
109 
110  for (const auto & old_node : mesh1.local_node_ptr_range())
111  {
112  const Node * new_node = *new_nit++;
113 
114  // Mesh::operator= hopefully preserved elem/node orderings...
115  libmesh_assert (*old_node == *new_node);
116 
117  for (int argi = 4; argi < argc; argi += 2)
118  {
119  const unsigned int pairnum = (argi-4)/2;
120 
121  const System & old_sys = *old_system[pairnum];
122  System & new_sys = *new_system[pairnum];
123 
124  const unsigned int n_comp =
125  old_node->n_comp(old_sys_num[pairnum],old_var_num[pairnum]);
126  libmesh_assert_equal_to(n_comp,
127  new_node->n_comp(new_sys_num[pairnum],new_var_num[pairnum]));
128 
129  for (unsigned int i=0; i<n_comp; i++)
130  {
131  const dof_id_type
132  old_index = old_node->dof_number
133  (old_sys_num[pairnum], old_var_num[pairnum], i),
134  new_index = new_node->dof_number
135  (new_sys_num[pairnum], new_var_num[pairnum], i);
136  new_sys.solution->set(new_index,(*old_sys.solution)(old_index));
137  }
138  }
139  }
140 
141 
142  // Copy over any element degree of freedom coefficients
144 
145  for (const auto & old_elem : mesh1.active_local_element_ptr_range())
146  {
147  const Elem * new_elem = *new_eit++;
148 
149  // Mesh::operator= hopefully preserved elem/node orderings...
150  libmesh_assert (*old_elem == *new_elem);
151 
152  for (int argi = 4; argi < argc; argi += 2)
153  {
154  const unsigned int pairnum = (argi-4)/2;
155 
156  const System & old_sys = *old_system[pairnum];
157  System & new_sys = *new_system[pairnum];
158 
159  const unsigned int n_comp =
160  old_elem->n_comp(old_sys_num[pairnum],old_var_num[pairnum]);
161  libmesh_assert_equal_to(n_comp,
162  new_elem->n_comp(new_sys_num[pairnum],new_var_num[pairnum]));
163 
164  for (unsigned int i=0; i<n_comp; i++)
165  {
166  const dof_id_type
167  old_index = old_elem->dof_number
168  (old_sys_num[pairnum], old_var_num[pairnum], i),
169  new_index = new_elem->dof_number
170  (new_sys_num[pairnum], new_var_num[pairnum], i);
171  new_sys.solution->set(new_index,(*old_sys.solution)(old_index));
172  }
173  }
174  }
175 
176  es2.write(argv[3], EquationSystems::WRITE_DATA);
177 
178  return 0;
179 }
int main(int argc, char **argv)
virtual element_iterator active_local_elements_begin() libmesh_override
const FEType & type() const
Definition: variable.h:119
This is the EquationSystems class.
A Node is like a Point, but with more information.
Definition: node.h:52
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.
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
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.
unsigned char dim
libmesh_assert(j)
This class defines the notion of a variable in the system.
Definition: variable.h:49
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.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
unsigned int number() const
Definition: system.h:2006
OStreamProxy out
The definition of the const_node_iterator struct.
Definition: mesh_base.h:1548
const std::set< subdomain_id_type > & active_subdomains() const
Definition: variable.h:150
virtual void init()
Initialize all the systems.
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:780
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual node_iterator local_nodes_begin() libmesh_override
Iterate over local nodes (nodes whose processor_id() matches the current processor).
uint8_t dof_id_type
Definition: id_types.h:64