Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
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 
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->num_cells = 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 =
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->num_cells = workset_size;
143 
144  if (!last_set_is_full) {
145  worksets.back().num_cells = last_workset_size;
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  Kokkos::View<int*,PHX::Device> cell_local_ids_k = Kokkos::View<int*,PHX::Device>("Workset:cell_local_ids",wkst->cell_local_ids.size());
158  for(std::size_t i=0;i<wkst->cell_local_ids.size();i++)
159  cell_local_ids_k(i) = wkst->cell_local_ids[i];
160  wkst->cell_local_ids_k = cell_local_ids_k;
161 
162  wkst->cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
163  vertex_coordinates.extent(1),
164  vertex_coordinates.extent(2));
165  wkst->block_id = elementBlock;
166  wkst->subcell_dim = needs.cellData.baseCellDimension();
167  wkst->subcell_index = 0;
168  }
169 
170  TEUCHOS_ASSERT(local_begin == local_cell_ids.end());
171 
172  // Copy cell vertex coordinates into local workset arrays
173  std::size_t offset = 0;
174  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
175  for (index_t cell = 0; cell < wkst->num_cells; ++cell)
176  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates.extent(1)); ++ vertex)
177  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates.extent(2)); ++ dim) {
178  //wkst->cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
179  wkst->cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
180  }
181 
182  offset += wkst->num_cells;
183  }
184 
185  TEUCHOS_ASSERT(offset == Teuchos::as<std::size_t>(vertex_coordinates.extent(0)));
186 
187  // Set ir and basis arrayskset
188  RCP<vector<int> > ir_degrees = rcp(new vector<int>(0));
189  RCP<vector<string> > basis_names = rcp(new vector<string>(0));
190  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
191  wkst->ir_degrees = ir_degrees;
192  wkst->basis_names = basis_names;
193  }
194 
195  // setup the integration rules and bases
196  for(std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst)
197  populateValueArrays(wkst->num_cells,false,needs,*wkst);
198 
199  return worksets_ptr;
200 }
201 
202 // ****************************************************************
203 // ****************************************************************
204 
205 template<typename ArrayT>
207 panzer::buildBCWorkset(const WorksetNeeds & needs,
208  const std::string & elementBlock,
209  const std::vector<std::size_t>& local_cell_ids,
210  const std::vector<std::size_t>& local_side_ids,
211  const ArrayT& vertex_coordinates,
212  const bool populate_value_arrays)
213 {
214  using Teuchos::RCP;
215  using Teuchos::rcp;
216 
217  panzer::MDFieldArrayFactory mdArrayFactory("",true);
218 
219  // key is local face index, value is workset with all elements
220  // for that local face
222  Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
223 
224  // All elements of boundary condition should go into one workset.
225  // However due to design of Intrepid2 (requires same basis for all
226  // cells), we have to separate the workset based on the local side
227  // index. Each workset for a boundary condition is associated with
228  // a local side for the element
229 
230  TEUCHOS_ASSERT(local_side_ids.size() == local_cell_ids.size());
231  TEUCHOS_ASSERT(local_side_ids.size() == static_cast<std::size_t>(vertex_coordinates.extent(0)));
232 
233  // key is local face index, value is a pair of cell index and vector of element local ids
234  std::map<unsigned,std::vector<std::pair<std::size_t,std::size_t> > > element_list;
235  for (std::size_t cell=0; cell < local_cell_ids.size(); ++cell)
236  element_list[local_side_ids[cell]].push_back(std::make_pair(cell,local_cell_ids[cell]));
237 
238  std::map<unsigned,panzer::Workset>& worksets = *worksets_ptr;
239 
240  // create worksets
241  std::map<unsigned,std::vector<std::pair<std::size_t,std::size_t> > >::const_iterator side;
242  for (side = element_list.begin(); side != element_list.end(); ++side) {
243 
244  std::vector<std::size_t>& cell_local_ids = worksets[side->first].cell_local_ids;
245 
246  worksets[side->first].cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",
247  side->second.size(),
248  vertex_coordinates.extent(1),
249  vertex_coordinates.extent(2));
250  Workset::CellCoordArray coords = worksets[side->first].cell_vertex_coordinates;
251 
252  for (std::size_t cell = 0; cell < side->second.size(); ++cell) {
253  cell_local_ids.push_back(side->second[cell].second);
254 
255  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates.extent(1)); ++ vertex)
256  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates.extent(2)); ++ dim) {
257  coords(cell,vertex,dim) = vertex_coordinates(side->second[cell].first,vertex,dim);
258  }
259  }
260 
261  Kokkos::View<int*,PHX::Device> cell_local_ids_k = Kokkos::View<int*,PHX::Device>("Workset:cell_local_ids",worksets[side->first].cell_local_ids.size());
262  for(std::size_t i=0;i<worksets[side->first].cell_local_ids.size();i++)
263  cell_local_ids_k(i) = worksets[side->first].cell_local_ids[i];
264  worksets[side->first].cell_local_ids_k = cell_local_ids_k;
265 
266  worksets[side->first].num_cells = worksets[side->first].cell_local_ids.size();
267  worksets[side->first].block_id = elementBlock;
268  worksets[side->first].subcell_dim = needs.cellData.baseCellDimension() - 1;
269  worksets[side->first].subcell_index = side->first;
270  }
271 
272  if (populate_value_arrays) {
273  // setup the integration rules and bases
274  for (std::map<unsigned,panzer::Workset>::iterator wkst = worksets.begin();
275  wkst != worksets.end(); ++wkst) {
276 
277  populateValueArrays(wkst->second.num_cells,true,needs,wkst->second); // populate "side" values
278  }
279  }
280 
281  return worksets_ptr;
282 }
283 
284 // ****************************************************************
285 // ****************************************************************
286 
287 namespace panzer {
288 namespace impl {
289 
290 /* Associate two sets of local side IDs into lists. Each list L has the property
291  * that every local side id in that list is the same, and this holds for each
292  * local side ID set. The smallest set of lists is found. The motivation for
293  * this procedure is to find a 1-1 workset pairing in advance. See the comment
294  * re: Intrepid2 in buildBCWorkset for more.
295  * The return value is an RCP to a map. Only the map's values are of interest
296  * in practice. Each value is a list L. The map's key is a pair (side ID a, side
297  * ID b) that gives rise to the list. We return a pointer to a map so that the
298  * caller does not have to deal with the map type; auto suffices.
299  */
301 associateCellsBySideIds(const std::vector<std::size_t>& sia /* local_side_ids_a */,
302  const std::vector<std::size_t>& sib /* local_side_ids_b */)
303 {
304  TEUCHOS_ASSERT(sia.size() == sib.size());
305 
306  auto sip2i_p = Teuchos::rcp(new std::map< std::pair<std::size_t, std::size_t>, std::vector<std::size_t> >);
307  auto& sip2i = *sip2i_p;
308 
309  for (std::size_t i = 0; i < sia.size(); ++i)
310  sip2i[std::make_pair(sia[i], sib[i])].push_back(i);
311 
312  return sip2i_p;
313 }
314 
315 // Set s = a(idxs). No error checking.
316 template <typename T>
317 void subset(const std::vector<T>& a, const std::vector<std::size_t>& idxs, std::vector<T>& s)
318 {
319  s.resize(idxs.size());
320  for (std::size_t i = 0; i < idxs.size(); ++i)
321  s[i] = a[idxs[i]];
322 }
323 
324 template<typename ArrayT>
327  const std::string & blockid_a,
328  const std::vector<std::size_t>& local_cell_ids_a,
329  const std::vector<std::size_t>& local_side_ids_a,
330  const ArrayT& vertex_coordinates_a,
331  const panzer::WorksetNeeds & needs_b,
332  const std::string & blockid_b,
333  const std::vector<std::size_t>& local_cell_ids_b,
334  const std::vector<std::size_t>& local_side_ids_b,
335  const ArrayT& vertex_coordinates_b,
336  const WorksetNeeds& needs_b2)
337 {
338  TEUCHOS_ASSERT(local_cell_ids_a.size() == local_cell_ids_b.size());
339  // Get a and b workset maps separately, but don't populate b's arrays.
341  mwa = buildBCWorkset(needs_a,blockid_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a),
342  mwb = buildBCWorkset(needs_b2, blockid_b, local_cell_ids_b, local_side_ids_b,
343  vertex_coordinates_b, false /* populate_value_arrays */);
344  TEUCHOS_ASSERT(mwa->size() == 1 && mwb->size() == 1);
345  for (std::map<unsigned,panzer::Workset>::iterator ait = mwa->begin(), bit = mwb->begin();
346  ait != mwa->end(); ++ait, ++bit) {
347  TEUCHOS_ASSERT(Teuchos::as<std::size_t>(ait->second.num_cells) == local_cell_ids_a.size() &&
348  Teuchos::as<std::size_t>(bit->second.num_cells) == local_cell_ids_b.size());
349  panzer::Workset& wa = ait->second;
350  // Copy b's details(0) to a's details(1).
351  wa.other = Teuchos::rcp(new panzer::WorksetDetails(bit->second.details(0)));
352  // Populate details(1) arrays so that IP are in order corresponding to details(0).
353  populateValueArrays(wa.num_cells, true, needs_b, wa.details(1), Teuchos::rcpFromRef(wa.details(0)));
354  }
355  // Now mwa has everything we need.
356  return mwa;
357 }
358 
359 } // namespace impl
360 } // namespace panzer
361 
362 // ****************************************************************
363 // ****************************************************************
364 
365 template<typename ArrayT>
368  const std::string & blockid_a,
369  const std::vector<std::size_t>& local_cell_ids_a,
370  const std::vector<std::size_t>& local_side_ids_a,
371  const ArrayT& vertex_coordinates_a,
372  const panzer::WorksetNeeds & needs_b,
373  const std::string & blockid_b,
374  const std::vector<std::size_t>& local_cell_ids_b,
375  const std::vector<std::size_t>& local_side_ids_b,
376  const ArrayT& vertex_coordinates_b)
377 {
378  // Since Intrepid2 requires all side IDs in a workset to be the same (see
379  // Intrepid2 comment above), we break the element list into pieces such that
380  // each piece contains elements on each side of the interface L_a and L_b and
381  // all elemnets L_a have the same side ID, and the same for L_b.
382  auto side_id_associations = impl::associateCellsBySideIds(local_side_ids_a, local_side_ids_b);
383  if (side_id_associations->size() == 1) {
384  // Common case of one workset on each side; optimize for it.
385  return impl::buildBCWorksetForUniqueSideId(needs_a, blockid_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a,
386  needs_b, blockid_b, local_cell_ids_b, local_side_ids_b, vertex_coordinates_b,
387  needs_b);
388  } else {
389  // The interface has elements having a mix of side IDs, so deal with each
390  // pair in turn.
391  Teuchos::RCP<std::map<unsigned,panzer::Workset> > mwa = Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
392  std::vector<std::size_t> lci_a, lci_b, lsi_a, lsi_b;
393  panzer::MDFieldArrayFactory mdArrayFactory("", true);
394  const int d1 = Teuchos::as<std::size_t>(vertex_coordinates_a.extent(1)),
395  d2 = Teuchos::as<std::size_t>(vertex_coordinates_a.extent(2));
396  for (auto it : *side_id_associations) {
397  const auto& idxs = it.second;
398  impl::subset(local_cell_ids_a, idxs, lci_a);
399  impl::subset(local_side_ids_a, idxs, lsi_a);
400  impl::subset(local_cell_ids_b, idxs, lci_b);
401  impl::subset(local_side_ids_b, idxs, lsi_b);
402  auto vc_a = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("vc_a", idxs.size(), d1, d2);
403  auto vc_b = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("vc_b", idxs.size(), d1, d2);
404  for (std::size_t i = 0; i < idxs.size(); ++i) {
405  const auto ii = idxs[i];
406  for (int j = 0; j < d1; ++j)
407  for (int k = 0; k < d2; ++k) {
408  vc_a(i, j, k) = vertex_coordinates_a(ii, j, k);
409  vc_b(i, j, k) = vertex_coordinates_b(ii, j, k);
410  }
411  }
412  auto mwa_it = impl::buildBCWorksetForUniqueSideId(needs_a,blockid_a, lci_a, lsi_a, vc_a,
413  needs_b,blockid_b, lci_b, lsi_b, vc_b,
414  needs_b);
415  TEUCHOS_ASSERT(mwa_it->size() == 1);
416  // Form a unique key that encodes the pair (side ID a, side ID b). We
417  // abuse the key here in the sense that it is everywhere else understood
418  // to correspond to the side ID of the elements in the workset.
419  // 1000 is a number substantially larger than is needed for any element.
420  const std::size_t key = lsi_a[0] * 1000 + lsi_b[0];
421  (*mwa)[key] = mwa_it->begin()->second;
422  }
423  return mwa;
424  }
425 }
426 
427 // ****************************************************************
428 // ****************************************************************
429 
430 template<typename ArrayT>
432 panzer::buildEdgeWorksets(const WorksetNeeds & needs_a,
433  const std::string & eblock_a,
434  const std::vector<std::size_t>& local_cell_ids_a,
435  const std::vector<std::size_t>& local_side_ids_a,
436  const ArrayT& vertex_coordinates_a,
437  const WorksetNeeds & needs_b,
438  const std::string & eblock_b,
439  const std::vector<std::size_t>& local_cell_ids_b,
440  const std::vector<std::size_t>& local_side_ids_b,
441  const ArrayT& vertex_coordinates_b)
442 {
443  using Teuchos::RCP;
444  using Teuchos::rcp;
445 
446  panzer::MDFieldArrayFactory mdArrayFactory("",true);
447 
448  std::size_t total_num_cells_a = local_cell_ids_a.size();
449  std::size_t total_num_cells_b = local_cell_ids_b.size();
450 
451  TEUCHOS_ASSERT(total_num_cells_a==total_num_cells_b);
452  TEUCHOS_ASSERT(local_side_ids_a.size() == local_cell_ids_a.size());
453  TEUCHOS_ASSERT(local_side_ids_a.size() == static_cast<std::size_t>(vertex_coordinates_a.extent(0)));
454  TEUCHOS_ASSERT(local_side_ids_b.size() == local_cell_ids_b.size());
455  TEUCHOS_ASSERT(local_side_ids_b.size() == static_cast<std::size_t>(vertex_coordinates_b.extent(0)));
456 
457  std::size_t total_num_cells = total_num_cells_a;
458 
459  std::size_t workset_size = needs_a.cellData.numCells();
460 
462  Teuchos::rcp(new std::vector<panzer::Workset>);
463  std::vector<panzer::Workset>& worksets = *worksets_ptr;
464 
465  // special case for 0 elements!
466  if(total_num_cells==0) {
467 
468  // Setup integration rules and basis
469  RCP<std::vector<int> > ir_degrees = rcp(new std::vector<int>(0));
470  RCP<std::vector<std::string> > basis_names = rcp(new std::vector<std::string>(0));
471 
472  worksets.resize(1);
473  std::vector<panzer::Workset>::iterator i = worksets.begin();
474 
475  i->details(0).block_id = eblock_a;
476  i->other = Teuchos::rcp(new panzer::WorksetDetails);
477  i->details(1).block_id = eblock_b;
478 
479  i->num_cells = 0;
480  i->ir_degrees = ir_degrees;
481  i->basis_names = basis_names;
482 
483  for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
484 
485  RCP<panzer::IntegrationValues2<double> > iv2 =
487  iv2->setupArrays(needs_a.int_rules[j]);
488 
489  ir_degrees->push_back(needs_a.int_rules[j]->cubature_degree);
490  i->int_rules.push_back(iv2);
491  }
492 
493  // Need to create all combinations of basis/ir pairings
494  for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
495 
496  for(std::size_t b=0;b<needs_a.bases.size();b++) {
497 
498  RCP<panzer::BasisIRLayout> b_layout = rcp(new panzer::BasisIRLayout(needs_a.bases[b],*needs_a.int_rules[j]));
499 
500  RCP<panzer::BasisValues2<double> > bv2 =
501  rcp(new panzer::BasisValues2<double>("",true,true));
502  bv2->setupArrays(b_layout);
503  i->bases.push_back(bv2);
504 
505  basis_names->push_back(b_layout->name());
506  }
507 
508  }
509 
510  return worksets_ptr;
511  } // end special case
512 
513  // This collects all the elements that share the same sub cell pairs, this makes it easier to
514  // build the required worksets
515  // key is the pair of local face indices, value is a vector of cell indices that satisfy this pair
516  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > element_list;
517  for (std::size_t cell=0; cell < local_cell_ids_a.size(); ++cell)
518  element_list[std::make_pair(local_side_ids_a[cell],local_side_ids_b[cell])].push_back(cell);
519 
520  // this is the lone iterator that will be used to loop over the element edge list
521  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator edge;
522 
523  // figure out how many worksets will be needed, resize workset vector accordingly
524  std::size_t num_worksets = 0;
525  for(edge=element_list.begin(); edge!=element_list.end();++edge) {
526  std::size_t num_worksets_for_edge = edge->second.size() / workset_size;
527  std::size_t last_workset_size = edge->second.size() % workset_size;
528  if(last_workset_size!=0)
529  num_worksets_for_edge += 1;
530 
531  num_worksets += num_worksets_for_edge;
532  }
533  worksets.resize(num_worksets);
534 
535  // fill the worksets
536  std::vector<Workset>::iterator current_workset = worksets.begin();
537  for(edge=element_list.begin(); edge!=element_list.end();++edge) {
538  // loop over each workset
539  const std::vector<std::size_t> & cell_indices = edge->second;
540 
541  current_workset = buildEdgeWorksets(cell_indices,
542  needs_a,eblock_a,local_cell_ids_a,local_side_ids_a,vertex_coordinates_a,
543  needs_b,eblock_b,local_cell_ids_b,local_side_ids_b,vertex_coordinates_b,
544  current_workset);
545  }
546 
547  // sanity check
548  TEUCHOS_ASSERT(current_workset==worksets.end());
549 
550  return worksets_ptr;
551 }
552 
553 template<typename ArrayT>
554 std::vector<panzer::Workset>::iterator
555 panzer::buildEdgeWorksets(const std::vector<std::size_t> & cell_indices,
556  const WorksetNeeds & needs_a,
557  const std::string & eblock_a,
558  const std::vector<std::size_t>& local_cell_ids_a,
559  const std::vector<std::size_t>& local_side_ids_a,
560  const ArrayT& vertex_coordinates_a,
561  const WorksetNeeds & needs_b,
562  const std::string & eblock_b,
563  const std::vector<std::size_t>& local_cell_ids_b,
564  const std::vector<std::size_t>& local_side_ids_b,
565  const ArrayT& vertex_coordinates_b,
566  std::vector<Workset>::iterator beg)
567 {
568  panzer::MDFieldArrayFactory mdArrayFactory("",true);
569 
570  std::vector<Workset>::iterator wkst = beg;
571 
572  std::size_t current_cell_index = 0;
573  while (current_cell_index<cell_indices.size()) {
574  std::size_t workset_size = needs_a.cellData.numCells();
575 
576  // allocate workset details (details(0) is already associated with the
577  // workset object itself)
578  wkst->other = Teuchos::rcp(new panzer::WorksetDetails);
579 
580  wkst->subcell_dim = needs_a.cellData.baseCellDimension()-1;
581 
582  wkst->details(0).subcell_index = local_side_ids_a[cell_indices[current_cell_index]];
583  wkst->details(0).block_id = eblock_a;
584  wkst->details(0).cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
585  vertex_coordinates_a.extent(1),
586  vertex_coordinates_a.extent(2));
587 
588  wkst->details(1).subcell_index = local_side_ids_b[cell_indices[current_cell_index]];
589  wkst->details(1).block_id = eblock_b;
590  wkst->details(1).cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
591  vertex_coordinates_a.extent(1),
592  vertex_coordinates_a.extent(2));
593 
594  std::size_t remaining_cells = cell_indices.size()-current_cell_index;
595  if(remaining_cells<workset_size)
596  workset_size = remaining_cells;
597 
598  // this is the true number of cells in this workset
599  wkst->num_cells = workset_size;
600  wkst->details(0).cell_local_ids.resize(workset_size);
601  wkst->details(1).cell_local_ids.resize(workset_size);
602 
603  for(std::size_t cell=0;cell<workset_size; cell++,current_cell_index++) {
604 
605  wkst->details(0).cell_local_ids[cell] = local_cell_ids_a[cell_indices[current_cell_index]];
606  wkst->details(1).cell_local_ids[cell] = local_cell_ids_b[cell_indices[current_cell_index]];
607 
608  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates_a.extent(1)); ++ vertex) {
609  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates_a.extent(2)); ++ dim) {
610  wkst->details(0).cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates_a(cell_indices[current_cell_index],vertex,dim);
611  wkst->details(1).cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates_b(cell_indices[current_cell_index],vertex,dim);
612  }
613  }
614  }
615 
616  Kokkos::View<int*,PHX::Device> cell_local_ids_k_0 = Kokkos::View<int*,PHX::Device>("Workset:cell_local_ids",wkst->details(0).cell_local_ids.size());
617  Kokkos::View<int*,PHX::Device> cell_local_ids_k_1 = Kokkos::View<int*,PHX::Device>("Workset:cell_local_ids",wkst->details(1).cell_local_ids.size());
618  for(std::size_t i=0;i<wkst->details(0).cell_local_ids.size();i++)
619  cell_local_ids_k_0(i) = wkst->details(0).cell_local_ids[i];
620  for(std::size_t i=0;i<wkst->details(1).cell_local_ids.size();i++)
621  cell_local_ids_k_1(i) = wkst->details(1).cell_local_ids[i];
622  wkst->details(0).cell_local_ids_k = cell_local_ids_k_0;
623  wkst->details(1).cell_local_ids_k = cell_local_ids_k_1;
624 
625  // fill the BasisValues and IntegrationValues arrays
626  std::size_t max_workset_size = needs_a.cellData.numCells();
627  populateValueArrays(max_workset_size,true,needs_a,wkst->details(0)); // populate "side" values
628  populateValueArrays(max_workset_size,true,needs_b,wkst->details(1),Teuchos::rcpFromRef(wkst->details(0)));
629 
630  wkst++;
631  }
632 
633  return wkst;
634 }
635 
636 #endif
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)
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
BASIS NODE
Spatial Dimension Tag.
#define TEUCHOS_ASSERT(assertion_test)
void populateValueArrays(std::size_t num_cells, bool isSide, const WorksetNeeds &needs, WorksetDetails &details, const Teuchos::RCP< WorksetDetails > other_details)