Panzer  Version of the Day
Panzer_Workset_Builder_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 #ifndef __Panzer_Workset_Builder_impl_hpp__
44 #define __Panzer_Workset_Builder_impl_hpp__
45 
46 #include <iostream>
47 #include <vector>
48 #include <map>
49 #include <algorithm>
50 #include "Panzer_Workset.hpp"
51 #include "Panzer_CellData.hpp"
52 #include "Panzer_BC.hpp"
55 
56 #include "Phalanx_DataLayout_MDALayout.hpp"
57 
58 // Intrepid2
59 #include "Shards_CellTopology.hpp"
60 #include "Intrepid2_DefaultCubatureFactory.hpp"
61 #include "Intrepid2_CellTools.hpp"
62 #include "Intrepid2_FunctionSpaceTools.hpp"
63 #include "Intrepid2_Basis.hpp"
64 
65 template<typename ArrayT>
66 Teuchos::RCP< std::vector<panzer::Workset> >
67 panzer::buildWorksets(const WorksetNeeds & needs,
68  const std::string & elementBlock,
69  const std::vector<std::size_t>& local_cell_ids,
70  const ArrayT& vertex_coordinates)
71 {
72  using std::vector;
73  using std::string;
74  using Teuchos::RCP;
75  using Teuchos::rcp;
76 
77  panzer::MDFieldArrayFactory mdArrayFactory("",true);
78 
79  std::size_t total_num_cells = local_cell_ids.size();
80 
81  std::size_t workset_size = needs.cellData.numCells();
82 
83  Teuchos::RCP< std::vector<panzer::Workset> > worksets_ptr =
84  Teuchos::rcp(new std::vector<panzer::Workset>);
85  std::vector<panzer::Workset>& worksets = *worksets_ptr;
86 
87  // special case for 0 elements!
88  if(total_num_cells==0) {
89 
90  // Setup integration rules and basis
91  RCP<vector<int> > ir_degrees = rcp(new vector<int>(0));
92  RCP<vector<string> > basis_names = rcp(new vector<string>(0));
93 
94  worksets.resize(1);
95  std::vector<panzer::Workset>::iterator i = worksets.begin();
96  i->setNumberOfCells(0,0,0);
97  i->block_id = elementBlock;
98  i->ir_degrees = ir_degrees;
99  i->basis_names = basis_names;
100 
101  for (std::size_t j=0;j<needs.int_rules.size();j++) {
102 
103  RCP<panzer::IntegrationValues2<double> > iv2 =
104  rcp(new panzer::IntegrationValues2<double>("",true));
105  iv2->setupArrays(needs.int_rules[j]);
106 
107  ir_degrees->push_back(needs.int_rules[j]->cubature_degree);
108  i->int_rules.push_back(iv2);
109  }
110 
111  // Need to create all combinations of basis/ir pairings
112  for (std::size_t j=0;j<needs.int_rules.size();j++) {
113  for (std::size_t b=0;b<needs.bases.size();b++) {
114  RCP<panzer::BasisIRLayout> b_layout
115  = rcp(new panzer::BasisIRLayout(needs.bases[b],*needs.int_rules[j]));
116 
117  RCP<panzer::BasisValues2<double> > bv2
118  = rcp(new panzer::BasisValues2<double>("",true,true));
119  bv2->setupArrays(b_layout);
120  i->bases.push_back(bv2);
121 
122  basis_names->push_back(b_layout->name());
123  }
124 
125  }
126 
127  return worksets_ptr;
128  } // end special case
129 
130  {
131  std::size_t num_worksets = total_num_cells / workset_size;
132  bool last_set_is_full = true;
133  std::size_t last_workset_size = total_num_cells % workset_size;
134  if (last_workset_size != 0) {
135  num_worksets += 1;
136  last_set_is_full = false;
137  }
138 
139  worksets.resize(num_worksets);
140  std::vector<panzer::Workset>::iterator i;
141  for (i = worksets.begin(); i != worksets.end(); ++i)
142  i->setNumberOfCells(workset_size,0,0);
143 
144  if (!last_set_is_full) {
145  worksets.back().setNumberOfCells(last_workset_size,0,0);
146  }
147  }
148 
149  // assign workset cell local ids
150  std::vector<std::size_t>::const_iterator local_begin = local_cell_ids.begin();
151  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
152  std::vector<std::size_t>::const_iterator begin_iter = local_begin;
153  std::vector<std::size_t>::const_iterator end_iter = begin_iter + wkst->num_cells;
154  local_begin = end_iter;
155  wkst->cell_local_ids.assign(begin_iter,end_iter);
156 
157  PHX::View<int*> cell_local_ids_k = PHX::View<int*>("Workset:cell_local_ids",wkst->cell_local_ids.size());
158  auto cell_local_ids_k_h = Kokkos::create_mirror_view(cell_local_ids_k);
159  for(std::size_t i=0;i<wkst->cell_local_ids.size();i++)
160  cell_local_ids_k_h(i) = wkst->cell_local_ids[i];
161  Kokkos::deep_copy(cell_local_ids_k, cell_local_ids_k_h);
162  wkst->cell_local_ids_k = cell_local_ids_k;
163 
164  wkst->cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
165  vertex_coordinates.extent(1),
166  vertex_coordinates.extent(2));
167  wkst->block_id = elementBlock;
168  wkst->subcell_dim = needs.cellData.baseCellDimension();
169  wkst->subcell_index = 0;
170  }
171 
172  TEUCHOS_ASSERT(local_begin == local_cell_ids.end());
173 
174  // Copy cell vertex coordinates into local workset arrays
175  std::size_t offset = 0;
176  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
177  auto cell_vertex_coordinates = wkst->cell_vertex_coordinates.get_static_view();
178  Kokkos::parallel_for(wkst->num_cells, KOKKOS_LAMBDA (int cell) {
179  for (std::size_t vertex = 0; vertex < vertex_coordinates.extent(1); ++ vertex)
180  for (std::size_t dim = 0; dim < vertex_coordinates.extent(2); ++ dim) {
181  //wkst->cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
182  cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
183  }
184  });
185  Kokkos::fence();
186  offset += wkst->num_cells;
187  }
188 
189  TEUCHOS_ASSERT(offset == Teuchos::as<std::size_t>(vertex_coordinates.extent(0)));
190 
191  // Set ir and basis arrayskset
192  RCP<vector<int> > ir_degrees = rcp(new vector<int>(0));
193  RCP<vector<string> > basis_names = rcp(new vector<string>(0));
194  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
195  wkst->ir_degrees = ir_degrees;
196  wkst->basis_names = basis_names;
197  }
198 
199  // setup the integration rules and bases
200  for(std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst)
201  populateValueArrays(wkst->num_cells,false,needs,*wkst);
202 
203  return worksets_ptr;
204 }
205 
206 // ****************************************************************
207 // ****************************************************************
208 
209 template<typename ArrayT>
210 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
211 panzer::buildBCWorkset(const WorksetNeeds & needs,
212  const std::string & elementBlock,
213  const std::vector<std::size_t>& local_cell_ids,
214  const std::vector<std::size_t>& local_side_ids,
215  const ArrayT& vertex_coordinates,
216  const bool populate_value_arrays)
217 {
218  using Teuchos::RCP;
219  using Teuchos::rcp;
220 
221  panzer::MDFieldArrayFactory mdArrayFactory("",true);
222 
223  // key is local face index, value is workset with all elements
224  // for that local face
225  auto worksets_ptr = Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
226 
227  // All elements of boundary condition should go into one workset.
228  // However due to design of Intrepid2 (requires same basis for all
229  // cells), we have to separate the workset based on the local side
230  // index. Each workset for a boundary condition is associated with
231  // a local side for the element
232 
233  TEUCHOS_ASSERT(local_side_ids.size() == local_cell_ids.size());
234  TEUCHOS_ASSERT(local_side_ids.size() == static_cast<std::size_t>(vertex_coordinates.extent(0)));
235 
236  // key is local face index, value is a pair of cell index and vector of element local ids
237  std::map< unsigned, std::vector< std::pair< std::size_t, std::size_t>>> element_list;
238  for (std::size_t cell = 0; cell < local_cell_ids.size(); ++cell)
239  element_list[local_side_ids[cell]].push_back(std::make_pair(cell, local_cell_ids[cell]));
240 
241  auto& worksets = *worksets_ptr;
242 
243  // create worksets
244  for (const auto& side : element_list) {
245 
246  auto& cell_local_ids = worksets[side.first].cell_local_ids;
247 
248  worksets[side.first].cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",
249  side.second.size(),
250  vertex_coordinates.extent(1),
251  vertex_coordinates.extent(2));
252  auto coords_view = worksets[side.first].cell_vertex_coordinates.get_view();
253  auto coords_h = Kokkos::create_mirror_view(coords_view);
254 
255  auto vertex_coordinates_h = Kokkos::create_mirror_view(PHX::as_view(vertex_coordinates));
256  Kokkos::deep_copy(vertex_coordinates_h, PHX::as_view(vertex_coordinates));
257 
258  for (std::size_t cell = 0; cell < side.second.size(); ++cell) {
259  cell_local_ids.push_back(side.second[cell].second);
260  const auto dim0 = side.second[cell].first;
261 
262  for(std::size_t vertex = 0; vertex < vertex_coordinates.extent(1); ++vertex)
263  {
264  const auto extent = Teuchos::as<std::size_t>(vertex_coordinates.extent(2));
265 
266  for (std::size_t dim = 0; dim < extent; ++dim)
267  coords_h(cell, vertex, dim) = vertex_coordinates_h(dim0, vertex, dim);
268  }
269  }
270 
271  Kokkos::deep_copy(coords_view, coords_h);
272 
273  const auto cell_local_ids_size = worksets[side.first].cell_local_ids.size();
274  auto cell_local_ids_k = PHX::View<int*>("Workset:cell_local_ids", cell_local_ids_size);
275  auto cell_local_ids_k_h = Kokkos::create_mirror_view(cell_local_ids_k);
276 
277  for(std::size_t i = 0; i < cell_local_ids_size; ++i){
278  cell_local_ids_k_h(i) = worksets.at(side.first).cell_local_ids[i];
279  }
280 
281  Kokkos::deep_copy(cell_local_ids_k, cell_local_ids_k_h);
282 
283  worksets[side.first].cell_local_ids_k = cell_local_ids_k;
284  worksets[side.first].num_cells = worksets[side.first].cell_local_ids.size();
285  worksets[side.first].block_id = elementBlock;
286  worksets[side.first].subcell_dim = needs.cellData.baseCellDimension() - 1;
287  worksets[side.first].subcell_index = side.first;
288  }
289 
290  if (populate_value_arrays) {
291  // setup the integration rules and bases
292  for (std::map<unsigned,panzer::Workset>::iterator wkst = worksets.begin();
293  wkst != worksets.end(); ++wkst) {
294 
295  populateValueArrays(wkst->second.num_cells,true,needs,wkst->second); // populate "side" values
296  }
297  }
298 
299  return worksets_ptr;
300 }
301 
302 // ****************************************************************
303 // ****************************************************************
304 
305 namespace panzer {
306 namespace impl {
307 
308 /* Associate two sets of local side IDs into lists. Each list L has the property
309  * that every local side id in that list is the same, and this holds for each
310  * local side ID set. The smallest set of lists is found. The motivation for
311  * this procedure is to find a 1-1 workset pairing in advance. See the comment
312  * re: Intrepid2 in buildBCWorkset for more.
313  * The return value is an RCP to a map. Only the map's values are of interest
314  * in practice. Each value is a list L. The map's key is a pair (side ID a, side
315  * ID b) that gives rise to the list. We return a pointer to a map so that the
316  * caller does not have to deal with the map type; auto suffices.
317  */
318 Teuchos::RCP< std::map<std::pair<std::size_t, std::size_t>, std::vector<std::size_t> > >
319 associateCellsBySideIds(const std::vector<std::size_t>& sia /* local_side_ids_a */,
320  const std::vector<std::size_t>& sib /* local_side_ids_b */)
321 {
322  TEUCHOS_ASSERT(sia.size() == sib.size());
323 
324  auto sip2i_p = Teuchos::rcp(new std::map< std::pair<std::size_t, std::size_t>, std::vector<std::size_t> >);
325  auto& sip2i = *sip2i_p;
326 
327  for (std::size_t i = 0; i < sia.size(); ++i)
328  sip2i[std::make_pair(sia[i], sib[i])].push_back(i);
329 
330  return sip2i_p;
331 }
332 
333 // Set s = a(idxs). No error checking.
334 template <typename T>
335 void subset(const std::vector<T>& a, const std::vector<std::size_t>& idxs, std::vector<T>& s)
336 {
337  s.resize(idxs.size());
338  for (std::size_t i = 0; i < idxs.size(); ++i)
339  s[i] = a[idxs[i]];
340 }
341 
342 template<typename ArrayT>
343 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
345  const std::string & blockid_a,
346  const std::vector<std::size_t>& local_cell_ids_a,
347  const std::vector<std::size_t>& local_side_ids_a,
348  const ArrayT& vertex_coordinates_a,
349  const panzer::WorksetNeeds & needs_b,
350  const std::string & blockid_b,
351  const std::vector<std::size_t>& local_cell_ids_b,
352  const std::vector<std::size_t>& local_side_ids_b,
353  const ArrayT& vertex_coordinates_b,
354  const WorksetNeeds& needs_b2)
355 {
356  TEUCHOS_ASSERT(local_cell_ids_a.size() == local_cell_ids_b.size());
357  // Get a and b workset maps separately, but don't populate b's arrays.
358  const Teuchos::RCP<std::map<unsigned,panzer::Workset> >
359  mwa = buildBCWorkset(needs_a,blockid_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a),
360  mwb = buildBCWorkset(needs_b2, blockid_b, local_cell_ids_b, local_side_ids_b,
361  vertex_coordinates_b, false /* populate_value_arrays */);
362  TEUCHOS_ASSERT(mwa->size() == 1 && mwb->size() == 1);
363  for (std::map<unsigned,panzer::Workset>::iterator ait = mwa->begin(), bit = mwb->begin();
364  ait != mwa->end(); ++ait, ++bit) {
365  TEUCHOS_ASSERT(Teuchos::as<std::size_t>(ait->second.num_cells) == local_cell_ids_a.size() &&
366  Teuchos::as<std::size_t>(bit->second.num_cells) == local_cell_ids_b.size());
367  panzer::Workset& wa = ait->second;
368  // Copy b's details(0) to a's details(1).
369  wa.other = Teuchos::rcp(new panzer::WorksetDetails(bit->second.details(0)));
370  // Populate details(1) arrays so that IP are in order corresponding to details(0).
371  populateValueArrays(wa.num_cells, true, needs_b, wa.details(1), Teuchos::rcpFromRef(wa.details(0)));
372  }
373  // Now mwa has everything we need.
374  return mwa;
375 }
376 
377 } // namespace impl
378 } // namespace panzer
379 
380 // ****************************************************************
381 // ****************************************************************
382 
383 template<typename ArrayT>
384 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
386  const std::string & blockid_a,
387  const std::vector<std::size_t>& local_cell_ids_a,
388  const std::vector<std::size_t>& local_side_ids_a,
389  const ArrayT& vertex_coordinates_a,
390  const panzer::WorksetNeeds & needs_b,
391  const std::string & blockid_b,
392  const std::vector<std::size_t>& local_cell_ids_b,
393  const std::vector<std::size_t>& local_side_ids_b,
394  const ArrayT& vertex_coordinates_b)
395 {
396  // Since Intrepid2 requires all side IDs in a workset to be the same (see
397  // Intrepid2 comment above), we break the element list into pieces such that
398  // each piece contains elements on each side of the interface L_a and L_b and
399  // all elemnets L_a have the same side ID, and the same for L_b.
400  auto side_id_associations = impl::associateCellsBySideIds(local_side_ids_a, local_side_ids_b);
401  if (side_id_associations->size() == 1) {
402  // Common case of one workset on each side; optimize for it.
403  return impl::buildBCWorksetForUniqueSideId(needs_a, blockid_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a,
404  needs_b, blockid_b, local_cell_ids_b, local_side_ids_b, vertex_coordinates_b,
405  needs_b);
406  } else {
407  // The interface has elements having a mix of side IDs, so deal with each
408  // pair in turn.
409  Teuchos::RCP<std::map<unsigned,panzer::Workset> > mwa = Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
410  std::vector<std::size_t> lci_a, lci_b, lsi_a, lsi_b;
411  panzer::MDFieldArrayFactory mdArrayFactory("", true);
412  const int d1 = Teuchos::as<std::size_t>(vertex_coordinates_a.extent(1)),
413  d2 = Teuchos::as<std::size_t>(vertex_coordinates_a.extent(2));
414  for (auto it : *side_id_associations) {
415  const auto& idxs = it.second;
416  impl::subset(local_cell_ids_a, idxs, lci_a);
417  impl::subset(local_side_ids_a, idxs, lsi_a);
418  impl::subset(local_cell_ids_b, idxs, lci_b);
419  impl::subset(local_side_ids_b, idxs, lsi_b);
420  auto vc_a = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("vc_a", idxs.size(), d1, d2);
421  auto vc_b = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("vc_b", idxs.size(), d1, d2);
422  auto vc_a_h = Kokkos::create_mirror_view(vc_a.get_static_view());
423  auto vc_b_h = Kokkos::create_mirror_view(vc_b.get_static_view());
424  auto vertex_coordinates_a_h = Kokkos::create_mirror_view(PHX::as_view(vertex_coordinates_a));
425  auto vertex_coordinates_b_h = Kokkos::create_mirror_view(PHX::as_view(vertex_coordinates_b));
426  Kokkos::deep_copy(vertex_coordinates_a_h, PHX::as_view(vertex_coordinates_a));
427  Kokkos::deep_copy(vertex_coordinates_b_h, PHX::as_view(vertex_coordinates_b));
428  for (std::size_t i = 0; i < idxs.size(); ++i) {
429  const auto ii = idxs[i];
430  for (int j = 0; j < d1; ++j)
431  for (int k = 0; k < d2; ++k) {
432  vc_a_h(i, j, k) = vertex_coordinates_a_h(ii, j, k);
433  vc_b_h(i, j, k) = vertex_coordinates_b_h(ii, j, k);
434  }
435  }
436  Kokkos::deep_copy(vc_a.get_static_view(), vc_a_h);
437  Kokkos::deep_copy(vc_b.get_static_view(), vc_b_h);
438  auto mwa_it = impl::buildBCWorksetForUniqueSideId(needs_a,blockid_a, lci_a, lsi_a, vc_a,
439  needs_b,blockid_b, lci_b, lsi_b, vc_b,
440  needs_b);
441  TEUCHOS_ASSERT(mwa_it->size() == 1);
442  // Form a unique key that encodes the pair (side ID a, side ID b). We
443  // abuse the key here in the sense that it is everywhere else understood
444  // to correspond to the side ID of the elements in the workset.
445  // 1000 is a number substantially larger than is needed for any element.
446  const std::size_t key = lsi_a[0] * 1000 + lsi_b[0];
447  (*mwa)[key] = mwa_it->begin()->second;
448  }
449  return mwa;
450  }
451 }
452 
453 // ****************************************************************
454 // ****************************************************************
455 
456 template<typename ArrayT>
457 Teuchos::RCP<std::vector<panzer::Workset> >
458 panzer::buildEdgeWorksets(const WorksetNeeds & needs_a,
459  const std::string & eblock_a,
460  const std::vector<std::size_t>& local_cell_ids_a,
461  const std::vector<std::size_t>& local_side_ids_a,
462  const ArrayT& vertex_coordinates_a,
463  const WorksetNeeds & needs_b,
464  const std::string & eblock_b,
465  const std::vector<std::size_t>& local_cell_ids_b,
466  const std::vector<std::size_t>& local_side_ids_b,
467  const ArrayT& vertex_coordinates_b)
468 {
469  using Teuchos::RCP;
470  using Teuchos::rcp;
471 
472  panzer::MDFieldArrayFactory mdArrayFactory("",true);
473 
474  std::size_t total_num_cells_a = local_cell_ids_a.size();
475  std::size_t total_num_cells_b = local_cell_ids_b.size();
476 
477  TEUCHOS_ASSERT(total_num_cells_a==total_num_cells_b);
478  TEUCHOS_ASSERT(local_side_ids_a.size() == local_cell_ids_a.size());
479  TEUCHOS_ASSERT(local_side_ids_a.size() == static_cast<std::size_t>(vertex_coordinates_a.extent(0)));
480  TEUCHOS_ASSERT(local_side_ids_b.size() == local_cell_ids_b.size());
481  TEUCHOS_ASSERT(local_side_ids_b.size() == static_cast<std::size_t>(vertex_coordinates_b.extent(0)));
482 
483  std::size_t total_num_cells = total_num_cells_a;
484 
485  std::size_t workset_size = needs_a.cellData.numCells();
486 
487  Teuchos::RCP< std::vector<panzer::Workset> > worksets_ptr =
488  Teuchos::rcp(new std::vector<panzer::Workset>);
489  std::vector<panzer::Workset>& worksets = *worksets_ptr;
490 
491  // special case for 0 elements!
492  if(total_num_cells==0) {
493 
494  // Setup integration rules and basis
495  RCP<std::vector<int> > ir_degrees = rcp(new std::vector<int>(0));
496  RCP<std::vector<std::string> > basis_names = rcp(new std::vector<std::string>(0));
497 
498  worksets.resize(1);
499  std::vector<panzer::Workset>::iterator i = worksets.begin();
500 
501  i->details(0).block_id = eblock_a;
502  i->other = Teuchos::rcp(new panzer::WorksetDetails);
503  i->details(1).block_id = eblock_b;
504 
505  i->num_cells = 0;
506  i->ir_degrees = ir_degrees;
507  i->basis_names = basis_names;
508 
509  for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
510 
511  RCP<panzer::IntegrationValues2<double> > iv2 =
512  rcp(new panzer::IntegrationValues2<double>("",true));
513  iv2->setupArrays(needs_a.int_rules[j]);
514 
515  ir_degrees->push_back(needs_a.int_rules[j]->cubature_degree);
516  i->int_rules.push_back(iv2);
517  }
518 
519  // Need to create all combinations of basis/ir pairings
520  for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
521 
522  for(std::size_t b=0;b<needs_a.bases.size();b++) {
523 
524  RCP<panzer::BasisIRLayout> b_layout = rcp(new panzer::BasisIRLayout(needs_a.bases[b],*needs_a.int_rules[j]));
525 
526  RCP<panzer::BasisValues2<double> > bv2 =
527  rcp(new panzer::BasisValues2<double>("",true,true));
528  bv2->setupArrays(b_layout);
529  i->bases.push_back(bv2);
530 
531  basis_names->push_back(b_layout->name());
532  }
533 
534  }
535 
536  return worksets_ptr;
537  } // end special case
538 
539  // This collects all the elements that share the same sub cell pairs, this makes it easier to
540  // build the required worksets
541  // key is the pair of local face indices, value is a vector of cell indices that satisfy this pair
542  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > element_list;
543  for (std::size_t cell=0; cell < local_cell_ids_a.size(); ++cell)
544  element_list[std::make_pair(local_side_ids_a[cell],local_side_ids_b[cell])].push_back(cell);
545 
546  // figure out how many worksets will be needed, resize workset vector accordingly
547  std::size_t num_worksets = 0;
548  for(const auto& edge : element_list) {
549  std::size_t num_worksets_for_edge = edge.second.size() / workset_size;
550  std::size_t last_workset_size = edge.second.size() % workset_size;
551  if(last_workset_size!=0)
552  num_worksets_for_edge += 1;
553 
554  num_worksets += num_worksets_for_edge;
555  }
556  worksets.resize(num_worksets);
557 
558  // fill the worksets
559  std::vector<Workset>::iterator current_workset = worksets.begin();
560  for(const auto& edge : element_list) {
561  // loop over each workset
562  const std::vector<std::size_t> & cell_indices = edge.second;
563 
564  current_workset = buildEdgeWorksets(cell_indices,
565  needs_a,eblock_a,local_cell_ids_a,local_side_ids_a,vertex_coordinates_a,
566  needs_b,eblock_b,local_cell_ids_b,local_side_ids_b,vertex_coordinates_b,
567  current_workset);
568  }
569 
570  // sanity check
571  TEUCHOS_ASSERT(current_workset==worksets.end());
572 
573  return worksets_ptr;
574 }
575 
576 template<typename ArrayT>
577 std::vector<panzer::Workset>::iterator
578 panzer::buildEdgeWorksets(const std::vector<std::size_t> & cell_indices,
579  const WorksetNeeds & needs_a,
580  const std::string & eblock_a,
581  const std::vector<std::size_t>& local_cell_ids_a,
582  const std::vector<std::size_t>& local_side_ids_a,
583  const ArrayT& vertex_coordinates_a,
584  const WorksetNeeds & needs_b,
585  const std::string & eblock_b,
586  const std::vector<std::size_t>& local_cell_ids_b,
587  const std::vector<std::size_t>& local_side_ids_b,
588  const ArrayT& vertex_coordinates_b,
589  std::vector<Workset>::iterator beg)
590 {
591  panzer::MDFieldArrayFactory mdArrayFactory("",true);
592 
593  std::vector<Workset>::iterator wkst = beg;
594 
595  std::size_t current_cell_index = 0;
596  while (current_cell_index<cell_indices.size()) {
597  std::size_t workset_size = needs_a.cellData.numCells();
598 
599  // allocate workset details (details(0) is already associated with the
600  // workset object itself)
601  wkst->other = Teuchos::rcp(new panzer::WorksetDetails);
602 
603  wkst->subcell_dim = needs_a.cellData.baseCellDimension()-1;
604 
605  wkst->details(0).subcell_index = local_side_ids_a[cell_indices[current_cell_index]];
606  wkst->details(0).block_id = eblock_a;
607  wkst->details(0).cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
608  vertex_coordinates_a.extent(1),
609  vertex_coordinates_a.extent(2));
610 
611  wkst->details(1).subcell_index = local_side_ids_b[cell_indices[current_cell_index]];
612  wkst->details(1).block_id = eblock_b;
613  wkst->details(1).cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
614  vertex_coordinates_a.extent(1),
615  vertex_coordinates_a.extent(2));
616 
617  std::size_t remaining_cells = cell_indices.size()-current_cell_index;
618  if(remaining_cells<workset_size)
619  workset_size = remaining_cells;
620 
621  // this is the true number of cells in this workset
622  wkst->setNumberOfCells(workset_size,0,0);
623  wkst->details(0).cell_local_ids.resize(workset_size);
624  wkst->details(1).cell_local_ids.resize(workset_size);
625 
626  auto dim0_cell_vertex_coordinates_view = wkst->details(0).cell_vertex_coordinates.get_static_view();
627  auto dim0_cell_vertex_coordinates_h = Kokkos::create_mirror_view(dim0_cell_vertex_coordinates_view);
628  Kokkos::deep_copy(dim0_cell_vertex_coordinates_h, dim0_cell_vertex_coordinates_view);
629 
630  auto dim1_cell_vertex_coordinates_view = wkst->details(1).cell_vertex_coordinates.get_static_view();
631  auto dim1_cell_vertex_coordinates_h = Kokkos::create_mirror_view(dim1_cell_vertex_coordinates_view);
632  Kokkos::deep_copy(dim1_cell_vertex_coordinates_h, dim1_cell_vertex_coordinates_view);
633 
634  auto vertex_coordinates_a_h = Kokkos::create_mirror_view(vertex_coordinates_a);
635  Kokkos::deep_copy(vertex_coordinates_a_h, vertex_coordinates_a);
636 
637  auto vertex_coordinates_b_h = Kokkos::create_mirror_view(vertex_coordinates_b);
638  Kokkos::deep_copy(vertex_coordinates_b_h, vertex_coordinates_b);
639 
640  for(std::size_t cell=0;cell<workset_size; cell++,current_cell_index++) {
641 
642  wkst->details(0).cell_local_ids[cell] = local_cell_ids_a[cell_indices[current_cell_index]];
643  wkst->details(1).cell_local_ids[cell] = local_cell_ids_b[cell_indices[current_cell_index]];
644 
645  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates_a.extent(1)); ++ vertex) {
646  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates_a.extent(2)); ++ dim) {
647  dim0_cell_vertex_coordinates_h(cell,vertex,dim) = vertex_coordinates_a_h(cell_indices[current_cell_index],vertex,dim);
648  dim1_cell_vertex_coordinates_h(cell,vertex,dim) = vertex_coordinates_b_h(cell_indices[current_cell_index],vertex,dim);
649  }
650  }
651  }
652 
653  Kokkos::deep_copy(dim0_cell_vertex_coordinates_view, dim0_cell_vertex_coordinates_h);
654  Kokkos::deep_copy(dim1_cell_vertex_coordinates_view, dim1_cell_vertex_coordinates_h);
655 
656  auto cell_local_ids_k_0 = PHX::View<int*>("Workset:cell_local_ids",wkst->details(0).cell_local_ids.size());
657  auto cell_local_ids_k_0_h = Kokkos::create_mirror_view(cell_local_ids_k_0);
658 
659  auto cell_local_ids_k_1 = PHX::View<int*>("Workset:cell_local_ids",wkst->details(1).cell_local_ids.size());
660  auto cell_local_ids_k_1_h = Kokkos::create_mirror_view(cell_local_ids_k_1);
661 
662  for(std::size_t i=0;i<wkst->details(0).cell_local_ids.size();i++)
663  cell_local_ids_k_0_h(i) = wkst->details(0).cell_local_ids[i];
664  for(std::size_t i=0;i<wkst->details(1).cell_local_ids.size();i++)
665  cell_local_ids_k_1_h(i) = wkst->details(1).cell_local_ids[i];
666 
667  Kokkos::deep_copy(cell_local_ids_k_0, cell_local_ids_k_0_h);
668  Kokkos::deep_copy(cell_local_ids_k_1, cell_local_ids_k_1_h);
669 
670  wkst->details(0).cell_local_ids_k = cell_local_ids_k_0;
671  wkst->details(1).cell_local_ids_k = cell_local_ids_k_1;
672 
673  // fill the BasisValues and IntegrationValues arrays
674  std::size_t max_workset_size = needs_a.cellData.numCells();
675  populateValueArrays(max_workset_size,true,needs_a,wkst->details(0)); // populate "side" values
676  populateValueArrays(max_workset_size,true,needs_b,wkst->details(1),Teuchos::rcpFromRef(wkst->details(0)));
677 
678  wkst++;
679  }
680 
681  return wkst;
682 }
683 
684 #endif
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
Teuchos::RCP< std::vector< Workset > > buildEdgeWorksets(const WorksetNeeds &needs_a, const std::string &eblock_a, const std::vector< std::size_t > &local_cell_ids_a, const std::vector< std::size_t > &local_side_ids_a, const ArrayT &vertex_coordinates_a, const WorksetNeeds &needs_b, const std::string &eblock_b, const std::vector< std::size_t > &local_cell_ids_b, const std::vector< std::size_t > &local_side_ids_b, const ArrayT &vertex_coordinates_b)
Teuchos::RCP< std::vector< Workset > > buildWorksets(const WorksetNeeds &needs, const std::string &elementBlock, const std::vector< std::size_t > &local_cell_ids, const ArrayT &vertex_coordinates)
Teuchos::RCP< std::map< unsigned, Workset > > buildBCWorkset(const WorksetNeeds &needs, const std::string &elementBlock, const std::vector< std::size_t > &local_cell_ids, const std::vector< std::size_t > &local_side_ids, const ArrayT &vertex_coordinates, const bool populate_value_arrays=true)
Teuchos::RCP< std::map< std::pair< std::size_t, std::size_t >, std::vector< std::size_t > > > associateCellsBySideIds(const std::vector< std::size_t > &sia, const std::vector< std::size_t > &sib)
Teuchos::RCP< WorksetDetails > other
Teuchos::RCP< std::map< unsigned, panzer::Workset > > buildBCWorksetForUniqueSideId(const panzer::WorksetNeeds &needs_a, const std::string &blockid_a, const std::vector< std::size_t > &local_cell_ids_a, const std::vector< std::size_t > &local_side_ids_a, const ArrayT &vertex_coordinates_a, const panzer::WorksetNeeds &needs_b, const std::string &blockid_b, const std::vector< std::size_t > &local_cell_ids_b, const std::vector< std::size_t > &local_side_ids_b, const ArrayT &vertex_coordinates_b, const WorksetNeeds &needs_b2)
void subset(const std::vector< T > &a, const std::vector< std::size_t > &idxs, std::vector< T > &s)
void populateValueArrays(std::size_t num_cells, bool isSide, const WorksetNeeds &needs, WorksetDetails &details, const Teuchos::RCP< WorksetDetails > other_details)