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