www.mooseframework.org
DTKInterpolationAdapter.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "libmesh/libmesh_config.h"
16 
17 #ifdef LIBMESH_TRILINOS_HAVE_DTK
18 
19 // MOOSE includes
20 #include "Moose.h"
23 #include "Transfer.h"
24 
25 #include "libmesh/mesh.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/equation_systems.h"
29 
30 // DTK includes
31 #include "libmesh/ignore_warnings.h"
32 #include <DTK_MeshTypes.hpp>
33 #include "libmesh/restore_warnings.h"
34 
35 DTKInterpolationAdapter::DTKInterpolationAdapter(Teuchos::RCP<const Teuchos::MpiComm<int>> in_comm,
36  EquationSystems & in_es,
37  const Point & offset,
38  unsigned int from_dim)
39  : comm(in_comm), es(in_es), _offset(offset), mesh(in_es.get_mesh()), dim(mesh.mesh_dimension())
40 {
41  Moose::ScopedCommSwapper swapper(*comm->getRawMpiComm());
42 
43  std::set<GlobalOrdinal> semi_local_nodes;
44  get_semi_local_nodes(semi_local_nodes);
45 
46  num_local_nodes = semi_local_nodes.size();
47 
48  vertices.resize(num_local_nodes);
49  Teuchos::ArrayRCP<double> coordinates(num_local_nodes * dim);
50 
51  Teuchos::ArrayRCP<double> target_coordinates(num_local_nodes * from_dim);
52 
53  // Fill in the vertices and coordinates
54  {
55  GlobalOrdinal i = 0;
56 
57  for (const auto & dof : semi_local_nodes)
58  {
59  const Node & node = mesh.node_ref(dof);
60 
61  vertices[i] = node.id();
62 
63  for (GlobalOrdinal j = 0; j < dim; j++)
64  coordinates[(j * num_local_nodes) + i] = node(j) + offset(j);
65 
66  for (GlobalOrdinal j = 0; j < from_dim; j++)
67  target_coordinates[(j * num_local_nodes) + i] = node(j) + offset(j);
68 
69  i++;
70  }
71  }
72 
73  // Currently assuming all elements are the same!
74  DataTransferKit::DTK_ElementTopology element_topology = get_element_topology(mesh.elem_ptr(0));
75  GlobalOrdinal n_nodes_per_elem = mesh.elem_ptr(0)->n_nodes();
76 
77  GlobalOrdinal n_local_elem = mesh.n_local_elem();
78 
79  elements.resize(n_local_elem);
80  Teuchos::ArrayRCP<GlobalOrdinal> connectivity(n_nodes_per_elem * n_local_elem);
81 
82  Teuchos::ArrayRCP<double> elem_centroid_coordinates(n_local_elem * from_dim);
83 
84  // Fill in the elements and connectivity
85  {
86  GlobalOrdinal i = 0;
87 
88  MeshBase::const_element_iterator end = mesh.local_elements_end();
89  for (MeshBase::const_element_iterator it = mesh.local_elements_begin(); it != end; ++it)
90  {
91  const Elem & elem = *(*it);
92  elements[i] = elem.id();
93 
94  for (GlobalOrdinal j = 0; j < n_nodes_per_elem; j++)
95  connectivity[(j * n_local_elem) + i] = elem.node_id(j);
96 
97  {
98  Point centroid = elem.centroid();
99  for (GlobalOrdinal j = 0; j < from_dim; j++)
100  elem_centroid_coordinates[(j * n_local_elem) + i] = centroid(j) + offset(j);
101  }
102 
103  i++;
104  }
105  }
106 
107  Teuchos::ArrayRCP<int> permutation_list(n_nodes_per_elem);
108  for (GlobalOrdinal i = 0; i < n_nodes_per_elem; ++i)
109  permutation_list[i] = i;
110 
111  /*
112  Moose::out<<"n_nodes_per_elem: "<<n_nodes_per_elem<<std::endl;
113 
114  Moose::out<<"Dim: "<<dim<<std::endl;
115 
116  Moose::err<<"Vertices size: "<<vertices.size()<<std::endl;
117  {
118  Moose::err<<libMesh::processor_id()<<" Vertices: ";
119 
120  for (unsigned int i=0; i<vertices.size(); i++)
121  Moose::err<<vertices[i]<<" ";
122 
123  Moose::err<<std::endl;
124  }
125 
126  Moose::err<<"Coordinates size: "<<coordinates.size()<<std::endl;
127  {
128  Moose::err<<libMesh::processor_id()<<" Coordinates: ";
129 
130  for (unsigned int i=0; i<coordinates.size(); i++)
131  Moose::err<<coordinates[i]<<" ";
132 
133  Moose::err<<std::endl;
134  }
135 
136  Moose::err<<"Connectivity size: "<<connectivity.size()<<std::endl;
137  {
138  Moose::err<<libMesh::processor_id()<<" Connectivity: ";
139 
140  for (unsigned int i=0; i<connectivity.size(); i++)
141  Moose::err<<connectivity[i]<<" ";
142 
143  Moose::err<<std::endl;
144  }
145 
146  Moose::err<<"Permutation_List size: "<<permutation_list.size()<<std::endl;
147  {
148  Moose::err<<libMesh::processor_id()<<" Permutation_List: ";
149 
150  for (unsigned int i=0; i<permutation_list.size(); i++)
151  Moose::err<<permutation_list[i]<<" ";
152 
153  Moose::err<<std::endl;
154  }
155 
156  */
157  Teuchos::RCP<MeshContainerType> mesh_container =
158  Teuchos::rcp(new MeshContainerType(dim,
159  vertices,
160  coordinates,
161  element_topology,
162  n_nodes_per_elem,
163  elements,
164  connectivity,
165  permutation_list));
166 
167  // We only have 1 element topology in this grid so we make just one mesh block
168  Teuchos::ArrayRCP<Teuchos::RCP<MeshContainerType>> mesh_blocks(1);
169  mesh_blocks[0] = mesh_container;
170 
171  // Create the MeshManager
172  mesh_manager =
173  Teuchos::rcp(new DataTransferKit::MeshManager<MeshContainerType>(mesh_blocks, comm, dim));
174 
175  // Pack the coordinates into a field, this will be the positions we'll ask for other systems
176  // fields at
177  if (from_dim == dim)
178  target_coords =
179  Teuchos::rcp(new DataTransferKit::FieldManager<MeshContainerType>(mesh_container, comm));
180  else
181  {
182  Teuchos::ArrayRCP<GlobalOrdinal> empty_elements(0);
183  Teuchos::ArrayRCP<GlobalOrdinal> empty_connectivity(0);
184 
185  Teuchos::RCP<MeshContainerType> coords_only_mesh_container =
186  Teuchos::rcp(new MeshContainerType(from_dim,
187  vertices,
188  target_coordinates,
189  element_topology,
190  n_nodes_per_elem,
191  empty_elements,
192  empty_connectivity,
193  permutation_list));
194 
195  target_coords = Teuchos::rcp(
196  new DataTransferKit::FieldManager<MeshContainerType>(coords_only_mesh_container, comm));
197  }
198 
199  {
200  Teuchos::ArrayRCP<GlobalOrdinal> empty_elements(0);
201  Teuchos::ArrayRCP<GlobalOrdinal> empty_connectivity(0);
202 
203  Teuchos::RCP<MeshContainerType> centroid_coords_only_mesh_container =
204  Teuchos::rcp(new MeshContainerType(from_dim,
205  elements,
206  elem_centroid_coordinates,
207  element_topology,
208  n_nodes_per_elem,
209  empty_elements,
210  empty_connectivity,
211  permutation_list));
212 
213  elem_centroid_coords = Teuchos::rcp(new DataTransferKit::FieldManager<MeshContainerType>(
214  centroid_coords_only_mesh_container, comm));
215  }
216 }
217 
220 {
221  if (evaluators.find(var_name) ==
222  evaluators.end()) // We haven't created an evaluator for the variable yet
223  {
224  System * sys = Transfer::find_sys(es, var_name);
225 
226  // Create the FieldEvaluator
227  evaluators[var_name] = Teuchos::rcp(new DTKInterpolationEvaluator(*sys, var_name, _offset));
228  }
229 
230  return evaluators[var_name];
231 }
232 
233 Teuchos::RCP<DataTransferKit::FieldManager<DTKInterpolationAdapter::FieldContainerType>>
235 {
236  if (values_to_fill.find(var_name) == values_to_fill.end())
237  {
238  System * sys = Transfer::find_sys(es, var_name);
239  unsigned int var_num = sys->variable_number(var_name);
240  bool is_nodal = sys->variable_type(var_num).family == LAGRANGE;
241 
242  Teuchos::ArrayRCP<double> data_space;
243 
244  if (is_nodal)
245  data_space = Teuchos::ArrayRCP<double>(vertices.size());
246  else
247  data_space = Teuchos::ArrayRCP<double>(elements.size());
248 
249  Teuchos::RCP<FieldContainerType> field_container =
250  Teuchos::rcp(new FieldContainerType(data_space, 1));
251  values_to_fill[var_name] =
252  Teuchos::rcp(new DataTransferKit::FieldManager<FieldContainerType>(field_container, comm));
253  }
254 
255  return values_to_fill[var_name];
256 }
257 
258 void
260  Teuchos::ArrayView<GlobalOrdinal> missed_points)
261 {
262  Moose::ScopedCommSwapper swapper(*comm->getRawMpiComm());
263 
264  System * sys = Transfer::find_sys(es, var_name);
265  unsigned int var_num = sys->variable_number(var_name);
266 
267  bool is_nodal = sys->variable_type(var_num).family == LAGRANGE;
268 
269  Teuchos::RCP<FieldContainerType> values = values_to_fill[var_name]->field();
270 
271  // Create a vector containing true or false for each point saying whether it was missed or not
272  // We're only going to update values for points that were not missed
273  std::vector<bool> missed(values->size(), false);
274 
275  for (const auto & dof : missed_points)
276  missed[dof] = true;
277 
278  unsigned int i = 0;
279  // Loop over the values (one for each node) and assign the value of this variable at each node
280  for (const auto & val : *values)
281  {
282  // If this point "missed" then skip it
283  if (missed[i])
284  {
285  i++;
286  continue;
287  }
288 
289  const DofObject * dof_object = NULL;
290 
291  if (is_nodal)
292  dof_object = mesh.node_ptr(vertices[i]);
293  else
294  dof_object = mesh.elem_ptr(elements[i]);
295 
296  if (dof_object->processor_id() == mesh.processor_id())
297  {
298  // The 0 is for the component... this only works for LAGRANGE!
299  dof_id_type dof = dof_object->dof_number(sys->number(), var_num, 0);
300  sys->solution->set(dof, val);
301  }
302 
303  i++;
304  }
305 
306  sys->solution->close();
307  sys->update();
308 }
309 
310 DataTransferKit::DTK_ElementTopology
312 {
313  ElemType type = elem->type();
314 
315  if (type == EDGE2)
316  return DataTransferKit::DTK_LINE_SEGMENT;
317  else if (type == EDGE3)
318  return DataTransferKit::DTK_LINE_SEGMENT;
319  else if (type == EDGE4)
320  return DataTransferKit::DTK_LINE_SEGMENT;
321  else if (type == TRI3)
322  return DataTransferKit::DTK_TRIANGLE;
323  else if (type == QUAD4)
324  return DataTransferKit::DTK_QUADRILATERAL;
325  else if (type == QUAD8)
326  return DataTransferKit::DTK_QUADRILATERAL;
327  else if (type == QUAD9)
328  return DataTransferKit::DTK_QUADRILATERAL;
329  else if (type == TET4)
330  return DataTransferKit::DTK_TETRAHEDRON;
331  else if (type == HEX8)
332  return DataTransferKit::DTK_HEXAHEDRON;
333  else if (type == HEX27)
334  return DataTransferKit::DTK_HEXAHEDRON;
335  else if (type == PYRAMID5)
336  return DataTransferKit::DTK_PYRAMID;
337 
338  libMesh::err << "Element type not supported by DTK!" << std::endl;
339  libmesh_error();
340 }
341 
342 void
343 DTKInterpolationAdapter::get_semi_local_nodes(std::set<GlobalOrdinal> & semi_local_nodes)
344 {
345  MeshBase::const_element_iterator end = mesh.local_elements_end();
346  for (MeshBase::const_element_iterator it = mesh.local_elements_begin(); it != end; ++it)
347  {
348  const Elem & elem = *(*it);
349 
350  for (unsigned int j = 0; j < elem.n_nodes(); j++)
351  semi_local_nodes.insert(elem.node_id(j));
352  }
353 }
354 
355 #endif // #ifdef LIBMESH_TRILINOS_HAVE_DTK
void update_variable_values(std::string var_name, Teuchos::ArrayView< GlobalOrdinal > missed_points)
After computing values for a variable in this EquationSystems we need to take those values and put th...
RCP_Evaluator get_variable_evaluator(std::string var_name)
std::map< std::string, Teuchos::RCP< DataTransferKit::FieldManager< FieldContainerType > > > values_to_fill
Map of variable names to arrays to be filled by a transfer.
DataTransferKit::FieldContainer< double > FieldContainerType
void get_semi_local_nodes(std::set< GlobalOrdinal > &semi_local_nodes)
Helper function that fills the std::set with all of the node numbers of nodes connected to local elem...
Teuchos::RCP< EvaluatorType > RCP_Evaluator
std::map< std::string, RCP_Evaluator > evaluators
Map of variable names to RCP_Evaluator objects.
DataTransferKit::MeshTraits< MeshContainerType >::global_ordinal_type GlobalOrdinal
Teuchos::ArrayRCP< GlobalOrdinal > vertices
Teuchos::RCP< DataTransferKit::MeshManager< MeshContainerType > > mesh_manager
DataTransferKit::DTK_ElementTopology get_element_topology(const libMesh::Elem *elem)
Helper that returns the DTK ElementTopology for a given Elem.
DataTransferKit::MeshContainer< long unsigned int > MeshContainerType
const libMesh::MeshBase & mesh
Teuchos::RCP< DataTransferKit::FieldManager< MeshContainerType > > elem_centroid_coords
Teuchos::RCP< DataTransferKit::FieldManager< MeshContainerType > > target_coords
Teuchos::RCP< const Teuchos::MpiComm< int > > comm
MatType type
libMesh::EquationSystems & es
MPI_Comm comm
Teuchos::ArrayRCP< GlobalOrdinal > elements
static System * find_sys(EquationSystems &es, const std::string &var_name)
Small helper function for finding the system containing the variable.
Definition: Transfer.C:70
Teuchos::RCP< DataTransferKit::FieldManager< FieldContainerType > > get_values_to_fill(std::string var_name)
DTKInterpolationAdapter(Teuchos::RCP< const Teuchos::MpiComm< int >> in_comm, libMesh::EquationSystems &in_es, const libMesh::Point &offset, unsigned int from_dim)