Panzer Version of the Day
Loading...
Searching...
No Matches
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
58
59#include <vector>
60#include <set>
61#include <string>
62
63namespace panzer
64{
65
66template <typename LocalOrdinal,typename GlobalOrdinal>
71
72#ifndef PANZER_HIDE_DEPRECATED_CODE
77template <typename LocalOrdinal,typename GlobalOrdinal>
83#endif
84
85template <typename LocalOrdinal,typename GlobalOrdinal>
88 const Teuchos::RCP<const Teuchos::Comm<int>> comm)
89{
90 initialize(conn, comm);
91}
92
93#ifndef PANZER_HIDE_DEPRECATED_CODE
98template <typename LocalOrdinal,typename GlobalOrdinal>
99void
102{
103 Teuchos::RCP<const Teuchos::Comm<int>> comm_world(new Teuchos::MpiComm< int>(MPI_COMM_WORLD));
104 initialize(conn, comm_world);
105}
106#endif
107
108template <typename LocalOrdinal,typename GlobalOrdinal>
109void
112 const Teuchos::RCP<const Teuchos::Comm<int>> comm)
113{
114 // Create a map of elems
115 std::vector<std::string> block_ids;
116 conn.getElementBlockIds(block_ids);
117
118 const int shift = 1;
119
120 int dimension;
121 std::vector<shards::CellTopology> ebt;
123 dimension = ebt[0].getDimension();
124
125 int my_rank = comm->getRank();
126#ifndef NDEBUG
127 int nprocs = comm->getSize();
128#endif
129
130 std::vector<GlobalOrdinal> element_GIDS;
131 for (size_t iblk = 0 ; iblk < block_ids.size(); ++iblk) {
132 // The connectivity takes in a shards class, therefore, it has to be build block by block?
133 // This seems odd, but o.k, moving forward.
134 if ( dimension == 1 ) {
135 panzer::EdgeFieldPattern edge_pattern(ebt[iblk]);
136 conn.buildConnectivity(edge_pattern);
137 } else if ( dimension == 2 ){
138 panzer::FaceFieldPattern face_pattern(ebt[iblk]);
139 conn.buildConnectivity(face_pattern);
140 } else {
141 panzer::ElemFieldPattern elem_pattern(ebt[iblk]);
142 conn.buildConnectivity(elem_pattern);
143 }
144 //const std::vector<GlobalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
145 const std::vector<LocalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
146 for (size_t i=0; i<block_elems.size(); ++i) {
147 const auto * connectivity = conn.getConnectivity(block_elems[i]);
148 element_GIDS.push_back(*connectivity);
149 }
150 }
151 Teuchos::RCP<const Map> elem_map =
152 Teuchos::RCP<Map>( new Map(Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(), &element_GIDS[0], element_GIDS.size(), 0, comm ));
153
154 // Now we need to create the face owned and owned/shared maps.
155 Teuchos::RCP<const Map> face_map,owned_face_map;
156 {
157 std::vector<GlobalOrdinal> face_GIDS;
158 std::set<GlobalOrdinal> set_of_face_GIDS;
159 for (size_t iblk = 0 ; iblk < block_ids.size(); ++iblk) {
160 // The connectivity takes in a shards class, therefore, it has to be build block by block?
161 // This seems odd, but o.k, moving forward.
162 if ( dimension == 1 ) {
163 panzer::NodalFieldPattern edge_pattern(ebt[iblk]);
164 conn.buildConnectivity(edge_pattern);
165 } else if ( dimension == 2 ){
166 panzer::EdgeFieldPattern face_pattern(ebt[iblk]);
167 conn.buildConnectivity(face_pattern);
168 } else {
169 panzer::FaceFieldPattern elem_pattern(ebt[iblk]);
170 conn.buildConnectivity(elem_pattern);
171 }
172 //const std::vector<GlobalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
173 const std::vector<LocalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
174 for (size_t i=0; i<block_elems.size(); ++i) {
175 int n_conn = conn.getConnectivitySize(block_elems[i]);
176 const panzer::GlobalOrdinal * connectivity = conn.getConnectivity(block_elems[i]);
177 for (int iface=0; iface<n_conn; ++iface)
178 set_of_face_GIDS.insert(connectivity[iface]);
179 }
180 }
181 face_GIDS.insert(face_GIDS.begin(), set_of_face_GIDS.begin(), set_of_face_GIDS.end());
182
183 face_map = Teuchos::RCP<Map>( new Map(Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(), &face_GIDS[0], face_GIDS.size(), 0, comm ));
184 owned_face_map = Tpetra::createOneToOne(face_map);
185 }
186
187 // OK, now we have a map of faces, and owned faces
188 // We are going to do create a multi vector of size 8, block/elem/proc/lidx, block/elem/proc/lidx
189 // (note lidx is the local index associted with a face on each cell)
190 Teuchos::RCP<GOMultiVector> owned_face2elem_mv = Teuchos::RCP<GOMultiVector>(new GOMultiVector(owned_face_map, 8));
191 Teuchos::RCP<GOMultiVector> face2elem_mv = Teuchos::RCP<GOMultiVector>(new GOMultiVector(face_map, 8));
192 // set a flag of -1-shift to identify unmodified data
193 face2elem_mv->putScalar(-1-shift);
194 {
195 auto f2e = face2elem_mv->getLocalViewHost(Tpetra::Access::ReadWrite);
196 auto b1 = Kokkos::subview(f2e,Kokkos::ALL(),0);
197 auto e1 = Kokkos::subview(f2e,Kokkos::ALL(),1);
198 auto p1 = Kokkos::subview(f2e,Kokkos::ALL(),2);
199 auto l1 = Kokkos::subview(f2e,Kokkos::ALL(),3);
200 auto b2 = Kokkos::subview(f2e,Kokkos::ALL(),4);
201 auto e2 = Kokkos::subview(f2e,Kokkos::ALL(),5);
202 auto p2 = Kokkos::subview(f2e,Kokkos::ALL(),6);
203 auto l2 = Kokkos::subview(f2e,Kokkos::ALL(),7);
204 // Now loop once again over the blocks
205 GlobalOrdinal my_elem = 0;
206 for (size_t iblk = 0 ; iblk < block_ids.size(); ++iblk) {
207 // The connectivity takes in a shards class, therefore, it has to be build block by block?
208 // This seems odd, but o.k, moving forward.
209 if ( dimension == 1 ) {
210 panzer::NodalFieldPattern edge_pattern(ebt[iblk]);
211 conn.buildConnectivity(edge_pattern);
212 } else if ( dimension == 2 ){
213 panzer::EdgeFieldPattern face_pattern(ebt[iblk]);
214 conn.buildConnectivity(face_pattern);
215 } else {
216 panzer::FaceFieldPattern elem_pattern(ebt[iblk]);
217 conn.buildConnectivity(elem_pattern);
218 }
219 //const std::vector<GlobalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
220 const std::vector<LocalOrdinal> &block_elems = conn.getElementBlock(block_ids[iblk]);
221 for (size_t i=0; i<block_elems.size(); ++i) {
222 int n_conn = conn.getConnectivitySize(block_elems[i]);
223 const panzer::GlobalOrdinal * connectivity = conn.getConnectivity(block_elems[i]);
224 for (int iface=0; iface<n_conn; ++iface) {
225 LocalOrdinal f = face_map->getLocalElement(connectivity[iface]);
226
227 // Fun fact: this method breaks if we have cell 0 found in block 0 on process 0 for local index 0
228 // Fix = add 'shift' to everything
229 if (b1[f] < 0 ) {
230 b1[f] = iblk+shift;
231 e1[f] = elem_map->getGlobalElement(my_elem)+shift;
232 p1[f] = my_rank+shift;
233 l1[f] = iface+shift;
234 } else if (b2[f] < 0){
235 b2[f] = iblk+shift;
236 e2[f] = elem_map->getGlobalElement(my_elem)+shift;
237 p2[f] = my_rank+shift;
238 l2[f] = iface+shift;
239 } else
240 assert(false);
241 }
242 ++my_elem;
243 }
244 }
245 }
246
247 // So, now we can export our owned things to our owned one.
248 Import imp(owned_face_map, face_map);
249 Export exp(face_map, owned_face_map);
250 owned_face2elem_mv->doExport(*face2elem_mv, exp, Tpetra::ADD);
251
252 {
253 auto f2e = face2elem_mv->getLocalViewHost(Tpetra::Access::ReadWrite);
254 auto b1 = Kokkos::subview(f2e,Kokkos::ALL(),0);
255 auto e1 = Kokkos::subview(f2e,Kokkos::ALL(),1);
256 auto p1 = Kokkos::subview(f2e,Kokkos::ALL(),2);
257 auto l1 = Kokkos::subview(f2e,Kokkos::ALL(),3);
258 auto b2 = Kokkos::subview(f2e,Kokkos::ALL(),4);
259 auto e2 = Kokkos::subview(f2e,Kokkos::ALL(),5);
260 auto p2 = Kokkos::subview(f2e,Kokkos::ALL(),6);
261 auto l2 = Kokkos::subview(f2e,Kokkos::ALL(),7);
262
263 auto of2e = owned_face2elem_mv->getLocalViewHost(Tpetra::Access::ReadWrite);
264 auto ob1 = Kokkos::subview(of2e,Kokkos::ALL(),0);
265 auto oe1 = Kokkos::subview(of2e,Kokkos::ALL(),1);
266 auto op1 = Kokkos::subview(of2e,Kokkos::ALL(),2);
267 auto ol1 = Kokkos::subview(of2e,Kokkos::ALL(),3);
268 auto ob2 = Kokkos::subview(of2e,Kokkos::ALL(),4);
269 auto oe2 = Kokkos::subview(of2e,Kokkos::ALL(),5);
270 auto op2 = Kokkos::subview(of2e,Kokkos::ALL(),6);
271 auto ol2 = Kokkos::subview(of2e,Kokkos::ALL(),7);
272
273 // Since we added all of the arrays together, they're going to be broken
274 // We need to fix all of the broken faces
275 int num_boundary=0;
276 for (size_t i=0; i<ob1.size();++i){
277
278 // Make sure side 1 of face was set (either by this process or by multiple processes
279 assert(b1[i] >= shift);
280
281 LocalOrdinal shared_local_id = face_map->getLocalElement(owned_face_map->getGlobalElement(i));
282 // handle purely internal faces
283 if (ob1[i] == b1[shared_local_id] && ob2[i] == b2[shared_local_id] &&
284 oe1[i] == e1[shared_local_id] && oe2[i] == e2[shared_local_id]) {
285 if (ob2[i] < 0 )
286 num_boundary++;
287 continue;
288 }
289
290 // Handle shared nodes on a boundary, this shouldn't happen
291 if (ob1[i] < b1[shared_local_id] || oe1[i] < e1[shared_local_id]) {
292 assert(false);
293 }
294
295 if ( ob1[i] > b1[shared_local_id] || oe1[i] > e1[shared_local_id]) {
296 // This case both wrote to a face, we need to detangle
297 assert(ob2[i] < 0 && oe2[i] < 0);
298 ob2[i] = ob1[i] - b1[shared_local_id];
299 oe2[i] = oe1[i] - e1[shared_local_id];
300 op2[i] = op1[i] - p1[shared_local_id];
301 ol2[i] = ol1[i] - l1[shared_local_id];
302
303 ob1[i] = b1[shared_local_id];
304 oe1[i] = e1[shared_local_id];
305 op1[i] = p1[shared_local_id];
306 ol1[i] = l1[shared_local_id];
307
308 assert(op1[i] >=0 && op2[i] >= 0 && op1[i] < nprocs+shift && op2[i] < nprocs+shift);
309 }
310 }
311 }
312 face2elem_mv->doImport(*owned_face2elem_mv, imp, Tpetra::REPLACE);
313
314 // std::cout << "number ext boundaries "<<num_boundary << std::endl;
315
316 // Now cache the data
317 face_map_ = face_map;
318 LocalOrdinal nfaces = face_map_->getLocalNumElements();
319 elems_by_face_ = PHX::View<GlobalOrdinal *[2]>("FaceToElement::elems_by_face_", nfaces);
320 lidx_by_face_ = PHX::View<int *[2]>("FaceToElement::elems_by_face_", nfaces);
321 blocks_by_face_ = PHX::View<int *[2]> ("FaceToElement::blocks_by_face_", nfaces);
322 procs_by_face_ = PHX::View<int *[2]> ("FaceToElement::procs_by_face_", nfaces);
323 auto elems_by_face_h = Kokkos::create_mirror_view(elems_by_face_);
324 auto blocks_by_face_h = Kokkos::create_mirror_view(blocks_by_face_);
325 auto procs_by_face_h = Kokkos::create_mirror_view(procs_by_face_);
326 auto lidx_by_face_h = Kokkos::create_mirror_view(lidx_by_face_);
327
328 auto f2e = face2elem_mv->getLocalViewHost(Tpetra::Access::ReadWrite);
329 auto b1 = Kokkos::subview(f2e,Kokkos::ALL(),0);
330 auto e1 = Kokkos::subview(f2e,Kokkos::ALL(),1);
331 auto p1 = Kokkos::subview(f2e,Kokkos::ALL(),2);
332 auto l1 = Kokkos::subview(f2e,Kokkos::ALL(),3);
333 auto b2 = Kokkos::subview(f2e,Kokkos::ALL(),4);
334 auto e2 = Kokkos::subview(f2e,Kokkos::ALL(),5);
335 auto p2 = Kokkos::subview(f2e,Kokkos::ALL(),6);
336 auto l2 = Kokkos::subview(f2e,Kokkos::ALL(),7);
337
338 for (LocalOrdinal i=0; i< nfaces; ++i) {
339 elems_by_face_h (i,0) = e1[i]-shift;
340 elems_by_face_h (i,1) = e2[i]-shift;
341 blocks_by_face_h(i,0) = b1[i]-shift;
342 blocks_by_face_h(i,1) = b2[i]-shift;
343 procs_by_face_h (i,0) = p1[i]-shift;
344 procs_by_face_h (i,1) = p2[i]-shift;
345 lidx_by_face_h (i,0) = l1[i]-shift;
346 lidx_by_face_h (i,1) = l2[i]-shift;
347
348 if(elems_by_face_h (i,0) < 0){
349 elems_by_face_h (i,0) = -1;
350 }
351 if(elems_by_face_h (i,1) < 0){
352 elems_by_face_h (i,1) = -1;
353 }
354 }
355 Kokkos::deep_copy(elems_by_face_, elems_by_face_h);
356 Kokkos::deep_copy(blocks_by_face_, blocks_by_face_h);
357 Kokkos::deep_copy(procs_by_face_, procs_by_face_h);
358 Kokkos::deep_copy(lidx_by_face_, lidx_by_face_h);
359}
360
361}
362
363#endif /* __FaceToElementent_impl_hpp__ */
Pure virtual base class for supplying mesh connectivity information to the DOF Manager.
virtual const std::vector< LocalOrdinal > & getElementBlock(const std::string &blockID) const =0
virtual void buildConnectivity(const FieldPattern &fp)=0
virtual const GlobalOrdinal * getConnectivity(LocalOrdinal localElmtId) const =0
virtual void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
virtual void getElementBlockTopologies(std::vector< shards::CellTopology > &elementBlockTopologies) const =0
virtual LocalOrdinal getConnectivitySize(LocalOrdinal localElmtId) const =0
Tpetra::MultiVector< GlobalOrdinal, LocalOrdinal, GlobalOrdinal, NodeType > GOMultiVector
Tpetra::Export< LocalOrdinal, GlobalOrdinal, NodeType > Export
void initialize(panzer::ConnManager &conn)
Tpetra::Map< LocalOrdinal, GlobalOrdinal, NodeType > Map
Tpetra::Import< LocalOrdinal, GlobalOrdinal, NodeType > Import