Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_FaceToElement_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 /*
44  * FaceToElement.cpp
45  *
46  * Created on: Nov 15, 2016
47  * Author: mbetten
48  */
49 
50 #ifndef PANZER_FACE_TO_ELEMENT_IMPL_HPP
51 #define PANZER_FACE_TO_ELEMENT_IMPL_HPP
52 
53 #include "Panzer_FaceToElement.hpp"
58 
59 #include <vector>
60 #include <set>
61 #include <string>
62 
63 namespace panzer
64 {
65 
66 template <typename LocalOrdinal,typename GlobalOrdinal>
69 {
70 }
71 
72 template <typename LocalOrdinal,typename GlobalOrdinal>
75 {
76  initialize(conn);
77 }
78 
79 template <typename LocalOrdinal,typename GlobalOrdinal>
80 void
83 {
84  // Create a map of elems
85  std::vector<std::string> block_ids;
86  conn.getElementBlockIds(block_ids);
87 
88  const int shift = 1;
89 
90  int dimension;
91  std::vector<shards::CellTopology> ebt;
92  conn.getElementBlockTopologies(ebt);
93  dimension = ebt[0].getDimension();
94 
96  int my_rank = comm->getRank();
97 #ifndef NDEBUG
98  int nprocs = comm->getSize();
99 #endif
100 
101  std::vector<GlobalOrdinal> element_GIDS;
102  for (size_t iblk = 0 ; iblk < block_ids.size(); ++iblk) {
103  // The connectivity takes in a shards class, therefore, it has to be build block by block?
104  // This seems odd, but o.k, moving forward.
105  if ( dimension == 1 ) {
106  panzer::EdgeFieldPattern edge_pattern(ebt[iblk]);
107  conn.buildConnectivity(edge_pattern);
108  } else if ( dimension == 2 ){
109  panzer::FaceFieldPattern face_pattern(ebt[iblk]);
110  conn.buildConnectivity(face_pattern);
111  } else {
112  panzer::ElemFieldPattern elem_pattern(ebt[iblk]);
113  conn.buildConnectivity(elem_pattern);
114  }
115  //const std::vector<GlobalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
116  const std::vector<LocalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
117  for (size_t i=0; i<block_elems.size(); ++i) {
118  const GlobalOrdinal * connectivity = conn.getConnectivity(block_elems[i]);
119  element_GIDS.push_back(*connectivity);
120  }
121  }
122  Teuchos::RCP<const Map> elem_map =
123  Teuchos::RCP<Map>( new Map(Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(), &element_GIDS[0], element_GIDS.size(), 0, comm ));
124 
125  // Now we need to create the face owned and owned/shared maps.
126  Teuchos::RCP<const Map> face_map,owned_face_map;
127  {
128  std::vector<GlobalOrdinal> face_GIDS;
129  std::set<GlobalOrdinal> set_of_face_GIDS;
130  for (size_t iblk = 0 ; iblk < block_ids.size(); ++iblk) {
131  // The connectivity takes in a shards class, therefore, it has to be build block by block?
132  // This seems odd, but o.k, moving forward.
133  if ( dimension == 1 ) {
134  panzer::NodalFieldPattern edge_pattern(ebt[iblk]);
135  conn.buildConnectivity(edge_pattern);
136  } else if ( dimension == 2 ){
137  panzer::EdgeFieldPattern face_pattern(ebt[iblk]);
138  conn.buildConnectivity(face_pattern);
139  } else {
140  panzer::FaceFieldPattern elem_pattern(ebt[iblk]);
141  conn.buildConnectivity(elem_pattern);
142  }
143  //const std::vector<GlobalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
144  const std::vector<LocalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
145  for (size_t i=0; i<block_elems.size(); ++i) {
146  int n_conn = conn.getConnectivitySize(block_elems[i]);
147  const GlobalOrdinal * connectivity = conn.getConnectivity(block_elems[i]);
148  for (int iface=0; iface<n_conn; ++iface)
149  set_of_face_GIDS.insert(connectivity[iface]);
150  }
151  }
152  face_GIDS.insert(face_GIDS.begin(), set_of_face_GIDS.begin(), set_of_face_GIDS.end());
153 
154  face_map = Teuchos::RCP<Map>( new Map(Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(), &face_GIDS[0], face_GIDS.size(), 0, comm ));
155  owned_face_map = Tpetra::createOneToOne(face_map);
156  }
157 
158  // OK, now we have a map of faces, and owned faces
159  // We are going to do create a multi vector of size 8, block/elem/proc/lidx, block/elem/proc/lidx
160  // (note lidx is the local index associted with a face on each cell)
161  Teuchos::RCP<GOMultiVector> owned_face2elem_mv = Teuchos::RCP<GOMultiVector>(new GOMultiVector(owned_face_map, 8));
163  // set a flag of -1-shift to identify unmodified data
164  face2elem_mv->putScalar(-1-shift);
165  auto b1 = face2elem_mv->getDataNonConst(0);
166  auto e1 = face2elem_mv->getDataNonConst(1);
167  auto p1 = face2elem_mv->getDataNonConst(2);
168  auto l1 = face2elem_mv->getDataNonConst(3);
169  auto b2 = face2elem_mv->getDataNonConst(4);
170  auto e2 = face2elem_mv->getDataNonConst(5);
171  auto p2 = face2elem_mv->getDataNonConst(6);
172  auto l2 = face2elem_mv->getDataNonConst(7);
173  // Now loop once again over the blocks
174  GlobalOrdinal my_elem = 0;
175  for (size_t iblk = 0 ; iblk < block_ids.size(); ++iblk) {
176  // The connectivity takes in a shards class, therefore, it has to be build block by block?
177  // This seems odd, but o.k, moving forward.
178  if ( dimension == 1 ) {
179  panzer::NodalFieldPattern edge_pattern(ebt[iblk]);
180  conn.buildConnectivity(edge_pattern);
181  } else if ( dimension == 2 ){
182  panzer::EdgeFieldPattern face_pattern(ebt[iblk]);
183  conn.buildConnectivity(face_pattern);
184  } else {
185  panzer::FaceFieldPattern elem_pattern(ebt[iblk]);
186  conn.buildConnectivity(elem_pattern);
187  }
188  //const std::vector<GlobalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
189  const std::vector<LocalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
190  for (size_t i=0; i<block_elems.size(); ++i) {
191  int n_conn = conn.getConnectivitySize(block_elems[i]);
192  const GlobalOrdinal * connectivity = conn.getConnectivity(block_elems[i]);
193  for (int iface=0; iface<n_conn; ++iface) {
194  LocalOrdinal f = face_map->getLocalElement(connectivity[iface]);
195 
196  // Fun fact: this method breaks if we have cell 0 found in block 0 on process 0 for local index 0
197  // Fix = add 'shift' to everything
198  if (b1[f] < 0 ) {
199  b1[f] = iblk+shift;
200  e1[f] = elem_map->getGlobalElement(my_elem)+shift;
201  p1[f] = my_rank+shift;
202  l1[f] = iface+shift;
203  } else if (b2[f] < 0){
204  b2[f] = iblk+shift;
205  e2[f] = elem_map->getGlobalElement(my_elem)+shift;
206  p2[f] = my_rank+shift;
207  l2[f] = iface+shift;
208  } else
209  assert(false);
210  }
211  ++my_elem;
212  }
213  }
214 
215  // So, now we can export our owned things to our owned one.
216  Import imp(owned_face_map, face_map);
217  Export exp(face_map, owned_face_map);
218  owned_face2elem_mv->doExport(*face2elem_mv, exp, Tpetra::ADD);
219 
220  auto ob1 = owned_face2elem_mv->getDataNonConst(0);
221  auto oe1 = owned_face2elem_mv->getDataNonConst(1);
222  auto op1 = owned_face2elem_mv->getDataNonConst(2);
223  auto ol1 = owned_face2elem_mv->getDataNonConst(3);
224  auto ob2 = owned_face2elem_mv->getDataNonConst(4);
225  auto oe2 = owned_face2elem_mv->getDataNonConst(5);
226  auto op2 = owned_face2elem_mv->getDataNonConst(6);
227  auto ol2 = owned_face2elem_mv->getDataNonConst(7);
228 
229  // Since we added all of the arrays together, they're going to be broken
230  // We need to fix all of the broken faces
231  int num_boundary=0;
232  for (int i=0; i<ob1.size();++i){
233 
234  // Make sure side 1 of face was set (either by this process or by multiple processes
235  assert(b1[i] >= shift);
236 
237  LocalOrdinal shared_local_id = face_map->getLocalElement(owned_face_map->getGlobalElement(i));
238  // handle purely internal faces
239  if (ob1[i] == b1[shared_local_id] && ob2[i] == b2[shared_local_id] &&
240  oe1[i] == e1[shared_local_id] && oe2[i] == e2[shared_local_id]) {
241  if (ob2[i] < 0 )
242  num_boundary++;
243  continue;
244  }
245 
246  // Handle shared nodes on a boundary, this shouldn't happen
247  if (ob1[i] < b1[shared_local_id] || oe1[i] < e1[shared_local_id]) {
248  assert(false);
249  }
250 
251  if ( ob1[i] > b1[shared_local_id] || oe1[i] > e1[shared_local_id]) {
252  // This case both wrote to a face, we need to detangle
253  assert(ob2[i] < 0 && oe2[i] < 0);
254  ob2[i] = ob1[i] - b1[shared_local_id];
255  oe2[i] = oe1[i] - e1[shared_local_id];
256  op2[i] = op1[i] - p1[shared_local_id];
257  ol2[i] = ol1[i] - l1[shared_local_id];
258 
259  ob1[i] = b1[shared_local_id];
260  oe1[i] = e1[shared_local_id];
261  op1[i] = p1[shared_local_id];
262  ol1[i] = l1[shared_local_id];
263 
264  assert(op1[i] >=0 && op2[i] >= 0 && op1[i] < nprocs+shift && op2[i] < nprocs+shift);
265  }
266  }
267  face2elem_mv->doImport(*owned_face2elem_mv, imp, Tpetra::REPLACE);
268 
269  // std::cout << "number ext boundaries "<<num_boundary << std::endl;
270 
271  // Now cache the data
272  face_map_ = face_map;
273  LocalOrdinal nfaces = face_map_->getNodeNumElements();
274  elems_by_face_ = Kokkos::View<GlobalOrdinal *[2]>("FaceToElement::elems_by_face_", nfaces);
275  lidx_by_face_ = Kokkos::View<int *[2]>("FaceToElement::elems_by_face_", nfaces);
276  blocks_by_face_ = Kokkos::View<int *[2]> ("FaceToElement::blocks_by_face_", nfaces);
277  procs_by_face_ = Kokkos::View<int *[2]> ("FaceToElement::procs_by_face_", nfaces);
278 
279  // We have to subtract 'shift' because we added 'shift' earlier to shift things away from 0
280  for (LocalOrdinal i=0; i< nfaces; ++i) {
281  elems_by_face_ (i,0) = e1[i]-shift;
282  elems_by_face_ (i,1) = e2[i]-shift;
283  blocks_by_face_(i,0) = b1[i]-shift;
284  blocks_by_face_(i,1) = b2[i]-shift;
285  procs_by_face_ (i,0) = p1[i]-shift;
286  procs_by_face_ (i,1) = p2[i]-shift;
287  lidx_by_face_ (i,0) = l1[i]-shift;
288  lidx_by_face_ (i,1) = l2[i]-shift;
289 
290  if(elems_by_face_ (i,0) < 0){
291  elems_by_face_ (i,0) = -1;
292  }
293  if(elems_by_face_ (i,1) < 0){
294  elems_by_face_ (i,1) = -1;
295  }
296  }
297 }
298 
299 }
300 
301 #endif /* __FaceToElementent_impl_hpp__ */
virtual void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
Tpetra::Import< LocalOrdinal, GlobalOrdinal, NodeType > Import
virtual void getElementBlockTopologies(std::vector< shards::CellTopology > &elementBlockTopologies) const =0
Tpetra::Map< LocalOrdinal, GlobalOrdinal, NodeType > Map
virtual LocalOrdinal getConnectivitySize(LocalOrdinal localElmtId) const =0
virtual const GlobalOrdinal * getConnectivity(LocalOrdinal localElmtId) const =0
virtual const std::vector< LocalOrdinal > & getElementBlock(const std::string &blockID) const =0
Tpetra::Export< LocalOrdinal, GlobalOrdinal, NodeType > Export
Tpetra::MultiVector< GlobalOrdinal, LocalOrdinal, GlobalOrdinal, NodeType > GOMultiVector
virtual void buildConnectivity(const FieldPattern &fp)=0
void initialize(panzer::ConnManager< LocalOrdinal, GlobalOrdinal > &conn)