Loading [MathJax]/extensions/tex2jax.js
libMesh
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
dtk_adapter.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 
20 #include "libmesh/libmesh_config.h"
21 
22 #ifdef LIBMESH_TRILINOS_HAVE_DTK
23 
24 // libMesh includes
25 #include "libmesh/dtk_adapter.h"
26 #include "libmesh/dtk_evaluator.h"
27 #include "libmesh/int_range.h"
28 #include "libmesh/mesh.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/elem.h"
31 #include "libmesh/equation_systems.h"
32 
33 // Trilinos includes
34 #include "libmesh/ignore_warnings.h"
35 #include <DTK_MeshTypes.hpp>
36 #include <Teuchos_Comm.hpp>
37 #include "libmesh/restore_warnings.h"
38 
39 // C++ includes
40 #include <vector>
41 #include <numeric>
42 
43 namespace libMesh
44 {
45 
46 DTKAdapter::DTKAdapter(Teuchos::RCP<const Teuchos::Comm<int>> in_comm, EquationSystems & in_es):
47  comm(in_comm),
48  es(in_es),
49  mesh(in_es.get_mesh()),
50  dim(mesh.mesh_dimension())
51 {
52  std::set<unsigned int> semi_local_nodes;
53  get_semi_local_nodes(semi_local_nodes);
54 
55  num_local_nodes = semi_local_nodes.size();
56 
57  vertices.resize(num_local_nodes);
58  Teuchos::ArrayRCP<double> coordinates(num_local_nodes * dim);
59 
60  // Fill in the vertices and coordinates
61  {
62  unsigned int i = 0;
63 
64  for (const auto & id : semi_local_nodes)
65  {
66  const Node & node = mesh.node_ref(id);
67 
68  vertices[i] = node.id();
69 
70  for (unsigned int j=0; j<dim; j++)
71  coordinates[(j*num_local_nodes) + i] = node(j);
72 
73  i++;
74  }
75  }
76 
77  // Currently assuming all elements are the same!
78  DataTransferKit::DTK_ElementTopology element_topology = get_element_topology(mesh.elem_ptr(0));
79  unsigned int n_nodes_per_elem = mesh.elem_ptr(0)->n_nodes();
80 
81  unsigned int n_local_elem = mesh.n_local_elem();
82 
83  Teuchos::ArrayRCP<int> elements(n_local_elem);
84  Teuchos::ArrayRCP<int> connectivity(n_nodes_per_elem*n_local_elem);
85 
86  // Fill in the elements and connectivity
87  {
88  unsigned int i = 0;
89 
90  for (const auto & elem : as_range(mesh.local_elements_begin(),
91  mesh.local_elements_end()))
92  {
93  elements[i] = elem->id();
94 
95  for (unsigned int j=0; j<n_nodes_per_elem; j++)
96  connectivity[(j*n_local_elem)+i] = elem->node_id(j);
97 
98  i++;
99  }
100  }
101 
102  Teuchos::ArrayRCP<int> permutation_list(n_nodes_per_elem);
103  std::iota(permutation_list.begin(), permutation_list.end(), 0);
104 
105  /*
106  if (this->processor_id() == 1)
107  sleep(1);
108 
109  libMesh::out<<"n_nodes_per_elem: "<<n_nodes_per_elem<<std::endl;
110 
111  libMesh::out<<"Dim: "<<dim<<std::endl;
112 
113  libMesh::err<<"Vertices size: "<<vertices.size()<<std::endl;
114  {
115  libMesh::err<<this->processor_id()<<" Vertices: ";
116 
117  for (std::size_t i=0; i<vertices.size(); i++)
118  libMesh::err<<vertices[i]<<" ";
119 
120  libMesh::err<<std::endl;
121  }
122 
123  libMesh::err<<"Coordinates size: "<<coordinates.size()<<std::endl;
124  {
125  libMesh::err<<this->processor_id()<<" Coordinates: ";
126 
127  for (std::size_t i=0; i<coordinates.size(); i++)
128  libMesh::err<<coordinates[i]<<" ";
129 
130  libMesh::err<<std::endl;
131  }
132 
133  libMesh::err<<"Connectivity size: "<<connectivity.size()<<std::endl;
134  {
135  libMesh::err<<this->processor_id()<<" Connectivity: ";
136 
137  for (std::size_t i=0; i<connectivity.size(); i++)
138  libMesh::err<<connectivity[i]<<" ";
139 
140  libMesh::err<<std::endl;
141  }
142 
143  libMesh::err<<"Permutation_List size: "<<permutation_list.size()<<std::endl;
144  {
145  libMesh::err<<this->processor_id()<<" Permutation_List: ";
146 
147  for (std::size_t i=0; i<permutation_list.size(); i++)
148  libMesh::err<<permutation_list[i]<<" ";
149 
150  libMesh::err<<std::endl;
151  }
152 
153  */
154  Teuchos::RCP<MeshContainerType>
155  mesh_container = Teuchos::rcp(new MeshContainerType(dim,
156  vertices,
157  coordinates,
158  element_topology,
159  n_nodes_per_elem,
160  elements,
161  connectivity,
162  permutation_list));
163 
164  // We only have 1 element topology in this grid so we make just one mesh block
165  Teuchos::ArrayRCP<Teuchos::RCP<MeshContainerType>> mesh_blocks(1);
166  mesh_blocks[0] = mesh_container;
167 
168  // Create the MeshManager
169  mesh_manager = Teuchos::rcp(new DataTransferKit::MeshManager<MeshContainerType>(mesh_blocks, comm, dim) );
170 
171  // Pack the coordinates into a field, this will be the positions we'll ask for other systems fields at
172  target_coords = Teuchos::rcp(new DataTransferKit::FieldManager<MeshContainerType>(mesh_container, comm));
173 }
174 
177 {
178  // We try emplacing a nullptr into the "evaluators" map.
179  auto [it, emplaced] = evaluators.emplace(var_name, nullptr);
180 
181  // If the emplace succeeded, that means it was a new entry in the
182  // map, so we need to actually construct the object.
183  if (emplaced)
184  {
185  System * sys = find_sys(var_name);
186  it->second = Teuchos::rcp(new DTKEvaluator(*sys, var_name));
187  }
188 
189  return it->second;
190 }
191 
192 Teuchos::RCP<DataTransferKit::FieldManager<DTKAdapter::FieldContainerType>>
193 DTKAdapter::get_values_to_fill(std::string var_name)
194 {
195  // We try emplacing a nullptr into the "values_to_fill" map.
196  auto [it, emplaced] = values_to_fill.emplace(var_name, nullptr);
197 
198  // If the emplace succeeded, that means it was a new entry in the
199  // map, so we need to actually construct the object.
200  if (emplaced)
201  {
202  Teuchos::ArrayRCP<double> data_space(num_local_nodes);
203  Teuchos::RCP<FieldContainerType> field_container = Teuchos::rcp(new FieldContainerType(data_space, 1));
204  it->second = Teuchos::rcp(new DataTransferKit::FieldManager<FieldContainerType>(field_container, comm));
205  }
206 
207  return it->second;
208 }
209 
210 void
212 {
213  System * sys = find_sys(var_name);
214  unsigned int var_num = sys->variable_number(var_name);
215 
216  Teuchos::RCP<FieldContainerType> values = values_to_fill[var_name]->field();
217 
218  unsigned int i=0;
219  // Loop over the values (one for each node) and assign the value of this variable at each node
220  for (const auto & value : *values)
221  {
222  unsigned int node_num = vertices[i];
223  const Node & node = mesh.node_ref(node_num);
224 
225  if (node.processor_id() == sys->processor_id())
226  {
227  // The 0 is for the component... this only works for LAGRANGE!
228  dof_id_type dof = node.dof_number(sys->number(), var_num, 0);
229  sys->solution->set(dof, value);
230  }
231 
232  i++;
233  }
234 
235  sys->solution->close();
236 }
237 
238 
244 System *
245 DTKAdapter::find_sys(std::string var_name)
246 {
247  System * sys = nullptr;
248 
249  // Find the system this variable is from
250  for (auto i : make_range(es.n_systems()))
251  {
252  if (es.get_system(i).has_variable(var_name))
253  {
254  sys = &es.get_system(i);
255  break;
256  }
257  }
258 
259  libmesh_assert(sys);
260 
261  return sys;
262 }
263 
264 DataTransferKit::DTK_ElementTopology
266 {
267  ElemType type = elem->type();
268 
269  if (type == EDGE2)
270  return DataTransferKit::DTK_LINE_SEGMENT;
271  else if (type == TRI3)
272  return DataTransferKit::DTK_TRIANGLE;
273  else if (type == QUAD4)
274  return DataTransferKit::DTK_QUADRILATERAL;
275  else if (type == TET4)
276  return DataTransferKit::DTK_TETRAHEDRON;
277  else if (type == HEX8)
278  return DataTransferKit::DTK_HEXAHEDRON;
279  else if (type == PYRAMID5)
280  return DataTransferKit::DTK_PYRAMID;
281 
282  libmesh_error_msg("Element type not supported by DTK!");
283 }
284 
285 void
286 DTKAdapter::get_semi_local_nodes(std::set<unsigned int> & semi_local_nodes)
287 {
288  for (const auto & elem : as_range(mesh.local_elements_begin(), mesh.local_elements_end()))
289  for (const Node & node : elem->node_ref_range())
290  semi_local_nodes.insert(node.id());
291 }
292 
293 } // namespace libMesh
294 
295 #endif // #ifdef LIBMESH_TRILINOS_HAVE_DTK
ElemType
Defines an enum for geometric element types.
This is the EquationSystems class.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
std::map< std::string, RCP_Evaluator > evaluators
Map of variable names to RCP_Evaluator objects.
Definition: dtk_adapter.h:110
A Node is like a Point, but with more information.
Definition: node.h:52
unsigned int n_systems() const
Teuchos::RCP< const Teuchos::Comm< int > > comm
Definition: dtk_adapter.h:94
unsigned int dim
unsigned int dim
Definition: dtk_adapter.h:97
void update_variable_values(std::string var_name)
After computing values for a variable in this EquationSystems we need to take those values and put th...
Definition: dtk_adapter.C:211
unsigned int num_local_nodes
Definition: dtk_adapter.h:99
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
The libMesh namespace provides an interface to certain functionality in the library.
Implements the evaluate() function to compute FE solution values at points requested by DTK...
Definition: dtk_evaluator.h:60
const T_sys & get_system(std::string_view name) const
DataTransferKit::MeshContainer< int > MeshContainerType
Definition: dtk_adapter.h:57
const MeshBase & mesh
Definition: dtk_adapter.h:96
dof_id_type n_local_elem() const
Definition: mesh_base.h:548
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
Definition: utility.h:229
DataTransferKit::DTK_ElementTopology get_element_topology(const Elem *elem)
Definition: dtk_adapter.C:265
unsigned int variable_number(std::string_view var) const
Definition: system.C:1609
EquationSystems & es
Definition: dtk_adapter.h:95
unsigned int number() const
Definition: system.h:2350
RCP_Evaluator get_variable_evaluator(std::string var_name)
Definition: dtk_adapter.C:176
dof_id_type id() const
Definition: dof_object.h:828
virtual unsigned int n_nodes() const =0
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
libmesh_assert(ctx)
DTKAdapter(Teuchos::RCP< const Teuchos::Comm< int >> in_comm, EquationSystems &in_es)
Definition: dtk_adapter.C:46
Teuchos::RCP< EvaluatorType > RCP_Evaluator
Definition: dtk_adapter.h:61
virtual const Elem * elem_ptr(const dof_id_type i) const =0
static const bool value
Definition: xdr_io.C:54
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
std::map< std::string, Teuchos::RCP< DataTransferKit::FieldManager< FieldContainerType > > > values_to_fill
Map of variable names to arrays to be filled by a transfer.
Definition: dtk_adapter.h:107
Teuchos::RCP< DataTransferKit::FieldManager< MeshContainerType > > target_coords
Definition: dtk_adapter.h:104
void get_semi_local_nodes(std::set< unsigned int > &semi_local_nodes)
Helper function that fills the std::set with all of the node numbers of nodes connected to local elem...
Definition: dtk_adapter.C:286
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
processor_id_type processor_id() const
Teuchos::RCP< DataTransferKit::MeshManager< MeshContainerType > > mesh_manager
Definition: dtk_adapter.h:102
Teuchos::ArrayRCP< int > vertices
Definition: dtk_adapter.h:100
processor_id_type processor_id() const
Definition: dof_object.h:905
virtual ElemType type() const =0
DataTransferKit::FieldContainer< double > FieldContainerType
Definition: dtk_adapter.h:58
System * find_sys(std::string var_name)
Small helper function for finding the system containing the variable.
Definition: dtk_adapter.C:245
uint8_t dof_id_type
Definition: id_types.h:67
Teuchos::RCP< DataTransferKit::FieldManager< FieldContainerType > > get_values_to_fill(std::string var_name)
Definition: dtk_adapter.C:193