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"
53 #include "Panzer_PhysicsBlock.hpp"
56 
57 #include "Phalanx_DataLayout_MDALayout.hpp"
58 
59 // Intrepid2
60 #include "Shards_CellTopology.hpp"
61 #include "Intrepid2_DefaultCubatureFactory.hpp"
62 #include "Intrepid2_CellTools.hpp"
63 #include "Intrepid2_FunctionSpaceTools.hpp"
64 #include "Intrepid2_Basis.hpp"
65 
66 template<typename ArrayT>
67 Teuchos::RCP< std::vector<panzer::Workset> >
69  const std::vector<std::size_t>& local_cell_ids,
70  const ArrayT& vertex_coordinates)
71 {
72  using Teuchos::RCP;
73  using Teuchos::rcp;
74 
75  WorksetNeeds needs;
76  needs.cellData = pb.cellData();
77 
78  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules = pb.getIntegrationRules();
79  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules.begin();
80  ir_itr != int_rules.end(); ++ir_itr)
81  needs.int_rules.push_back(ir_itr->second);
82 
83  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases= pb.getBases();
84  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases.begin();
85  b_itr != bases.end(); ++b_itr)
86  needs.bases.push_back(b_itr->second);
87 
88  return buildWorksets(needs,pb.elementBlockID(),local_cell_ids,vertex_coordinates);
89 }
90 
91 template<typename ArrayT>
92 Teuchos::RCP< std::vector<panzer::Workset> >
93 panzer::buildWorksets(const WorksetNeeds & needs,
94  const std::string & elementBlock,
95  const std::vector<std::size_t>& local_cell_ids,
96  const ArrayT& vertex_coordinates)
97 {
98  using std::vector;
99  using std::string;
100  using Teuchos::RCP;
101  using Teuchos::rcp;
102 
103  panzer::MDFieldArrayFactory mdArrayFactory("",true);
104 
105  std::size_t total_num_cells = local_cell_ids.size();
106 
107  std::size_t workset_size = needs.cellData.numCells();
108 
109  Teuchos::RCP< std::vector<panzer::Workset> > worksets_ptr =
110  Teuchos::rcp(new std::vector<panzer::Workset>);
111  std::vector<panzer::Workset>& worksets = *worksets_ptr;
112 
113  // special case for 0 elements!
114  if(total_num_cells==0) {
115 
116  // Setup integration rules and basis
117  RCP<vector<int> > ir_degrees = rcp(new vector<int>(0));
118  RCP<vector<string> > basis_names = rcp(new vector<string>(0));
119 
120  worksets.resize(1);
121  std::vector<panzer::Workset>::iterator i = worksets.begin();
122  i->num_cells = 0;
123  i->block_id = elementBlock;
124  i->ir_degrees = ir_degrees;
125  i->basis_names = basis_names;
126 
127  for (std::size_t j=0;j<needs.int_rules.size();j++) {
128 
129  RCP<panzer::IntegrationValues2<double> > iv2 =
130  rcp(new panzer::IntegrationValues2<double>("",true));
131  iv2->setupArrays(needs.int_rules[j]);
132 
133  ir_degrees->push_back(needs.int_rules[j]->cubature_degree);
134  i->int_rules.push_back(iv2);
135  }
136 
137  // Need to create all combinations of basis/ir pairings
138  for (std::size_t j=0;j<needs.int_rules.size();j++) {
139  for (std::size_t b=0;b<needs.bases.size();b++) {
140  RCP<panzer::BasisIRLayout> b_layout
141  = rcp(new panzer::BasisIRLayout(needs.bases[b],*needs.int_rules[j]));
142 
143  RCP<panzer::BasisValues2<double> > bv2
144  = rcp(new panzer::BasisValues2<double>("",true,true));
145  bv2->setupArrays(b_layout);
146  i->bases.push_back(bv2);
147 
148  basis_names->push_back(b_layout->name());
149  }
150 
151  }
152 
153  return worksets_ptr;
154  } // end special case
155 
156  {
157  std::size_t num_worksets = total_num_cells / workset_size;
158  bool last_set_is_full = true;
159  std::size_t last_workset_size = total_num_cells % workset_size;
160  if (last_workset_size != 0) {
161  num_worksets += 1;
162  last_set_is_full = false;
163  }
164 
165  worksets.resize(num_worksets);
166  std::vector<panzer::Workset>::iterator i;
167  for (i = worksets.begin(); i != worksets.end(); ++i)
168  i->num_cells = workset_size;
169 
170  if (!last_set_is_full) {
171  worksets.back().num_cells = last_workset_size;
172  }
173  }
174 
175  // assign workset cell local ids
176  std::vector<std::size_t>::const_iterator local_begin = local_cell_ids.begin();
177  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
178  std::vector<std::size_t>::const_iterator begin_iter = local_begin;
179  std::vector<std::size_t>::const_iterator end_iter = begin_iter + wkst->num_cells;
180  local_begin = end_iter;
181  wkst->cell_local_ids.assign(begin_iter,end_iter);
182 
183  Kokkos::View<int*,PHX::Device> cell_local_ids_k = Kokkos::View<int*,PHX::Device>("Workset:cell_local_ids",wkst->cell_local_ids.size());
184  for(std::size_t i=0;i<wkst->cell_local_ids.size();i++)
185  cell_local_ids_k(i) = wkst->cell_local_ids[i];
186  wkst->cell_local_ids_k = cell_local_ids_k;
187 
188  wkst->cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
189  vertex_coordinates.dimension(1),
190  vertex_coordinates.dimension(2));
191  wkst->block_id = elementBlock;
192  wkst->subcell_dim = needs.cellData.baseCellDimension();
193  wkst->subcell_index = 0;
194  }
195 
196  TEUCHOS_ASSERT(local_begin == local_cell_ids.end());
197 
198  // Copy cell vertex coordinates into local workset arrays
199  std::size_t offset = 0;
200  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
201  for (index_t cell = 0; cell < wkst->num_cells; ++cell)
202  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates.dimension(1)); ++ vertex)
203  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates.dimension(2)); ++ dim) {
204  //wkst->cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
205  wkst->cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
206  }
207 
208  offset += wkst->num_cells;
209  }
210 
211  TEUCHOS_ASSERT(offset == Teuchos::as<std::size_t>(vertex_coordinates.dimension(0)));
212 
213  // Set ir and basis arrayskset
214  RCP<vector<int> > ir_degrees = rcp(new vector<int>(0));
215  RCP<vector<string> > basis_names = rcp(new vector<string>(0));
216  for (std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst) {
217  wkst->ir_degrees = ir_degrees;
218  wkst->basis_names = basis_names;
219  }
220 
221  // setup the integration rules and bases
222  for(std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst)
223  populateValueArrays(wkst->num_cells,false,needs,*wkst);
224 
225  return worksets_ptr;
226 }
227 
228 // ****************************************************************
229 // ****************************************************************
230 
231 template<typename ArrayT>
232 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
234  const std::vector<std::size_t>& local_cell_ids,
235  const std::vector<std::size_t>& local_side_ids,
236  const ArrayT& vertex_coordinates)
237 {
238  using Teuchos::RCP;
239  using Teuchos::rcp;
240 
241  WorksetNeeds needs;
242  needs.cellData = pb.cellData();
243 
244  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules = pb.getIntegrationRules();
245  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules.begin();
246  ir_itr != int_rules.end(); ++ir_itr)
247  needs.int_rules.push_back(ir_itr->second);
248 
249  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases= pb.getBases();
250  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases.begin();
251  b_itr != bases.end(); ++b_itr)
252  needs.bases.push_back(b_itr->second);
253 
254  return buildBCWorkset(needs,pb.elementBlockID(),local_cell_ids,local_side_ids,vertex_coordinates);
255 }
256 
257 template<typename ArrayT>
258 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
259 panzer::buildBCWorkset(const WorksetNeeds & needs,
260  const std::string & elementBlock,
261  const std::vector<std::size_t>& local_cell_ids,
262  const std::vector<std::size_t>& local_side_ids,
263  const ArrayT& vertex_coordinates,
264  const bool populate_value_arrays)
265 {
266  using Teuchos::RCP;
267  using Teuchos::rcp;
268 
269  panzer::MDFieldArrayFactory mdArrayFactory("",true);
270 
271  // key is local face index, value is workset with all elements
272  // for that local face
273  Teuchos::RCP<std::map<unsigned,panzer::Workset> > worksets_ptr =
274  Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
275 
276  // All elements of boundary condition should go into one workset.
277  // However due to design of Intrepid2 (requires same basis for all
278  // cells), we have to separate the workset based on the local side
279  // index. Each workset for a boundary condition is associated with
280  // a local side for the element
281 
282  TEUCHOS_ASSERT(local_side_ids.size() == local_cell_ids.size());
283  TEUCHOS_ASSERT(local_side_ids.size() == static_cast<std::size_t>(vertex_coordinates.dimension(0)));
284 
285  // key is local face index, value is a pair of cell index and vector of element local ids
286  std::map<unsigned,std::vector<std::pair<std::size_t,std::size_t> > > element_list;
287  for (std::size_t cell=0; cell < local_cell_ids.size(); ++cell)
288  element_list[local_side_ids[cell]].push_back(std::make_pair(cell,local_cell_ids[cell]));
289 
290  std::map<unsigned,panzer::Workset>& worksets = *worksets_ptr;
291 
292  // create worksets
293  std::map<unsigned,std::vector<std::pair<std::size_t,std::size_t> > >::const_iterator side;
294  for (side = element_list.begin(); side != element_list.end(); ++side) {
295 
296  std::vector<std::size_t>& cell_local_ids = worksets[side->first].cell_local_ids;
297 
298  worksets[side->first].cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",
299  side->second.size(),
300  vertex_coordinates.dimension(1),
301  vertex_coordinates.dimension(2));
302  Workset::CellCoordArray coords = worksets[side->first].cell_vertex_coordinates;
303 
304  for (std::size_t cell = 0; cell < side->second.size(); ++cell) {
305  cell_local_ids.push_back(side->second[cell].second);
306 
307  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates.dimension(1)); ++ vertex)
308  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates.dimension(2)); ++ dim) {
309  coords(cell,vertex,dim) = vertex_coordinates(side->second[cell].first,vertex,dim);
310  }
311  }
312 
313  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());
314  for(std::size_t i=0;i<worksets[side->first].cell_local_ids.size();i++)
315  cell_local_ids_k(i) = worksets[side->first].cell_local_ids[i];
316  worksets[side->first].cell_local_ids_k = cell_local_ids_k;
317 
318  worksets[side->first].num_cells = worksets[side->first].cell_local_ids.size();
319  worksets[side->first].block_id = elementBlock;
320  worksets[side->first].subcell_dim = needs.cellData.baseCellDimension() - 1;
321  worksets[side->first].subcell_index = side->first;
322  }
323 
324  if (populate_value_arrays) {
325  // setup the integration rules and bases
326  for (std::map<unsigned,panzer::Workset>::iterator wkst = worksets.begin();
327  wkst != worksets.end(); ++wkst) {
328 
329  populateValueArrays(wkst->second.num_cells,true,needs,wkst->second); // populate "side" values
330  }
331  }
332 
333  return worksets_ptr;
334 }
335 
336 // ****************************************************************
337 // ****************************************************************
338 
339 template<typename ArrayT>
340 Teuchos::RCP< std::vector<panzer::Workset> >
342  const std::vector<std::size_t>& local_cell_ids_a,
343  const std::vector<std::size_t>& local_side_ids_a,
344  const ArrayT& vertex_coordinates_a,
345  const panzer::PhysicsBlock & pb_b,
346  const std::vector<std::size_t>& local_cell_ids_b,
347  const std::vector<std::size_t>& local_side_ids_b,
348  const ArrayT& vertex_coordinates_b)
349 {
350  using Teuchos::RCP;
351  using Teuchos::rcp;
352 
353  WorksetNeeds needs_a;
354  needs_a.cellData = pb_a.cellData();
355 
356  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules_a = pb_a.getIntegrationRules();
357  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules_a.begin();
358  ir_itr != int_rules_a.end(); ++ir_itr)
359  needs_a.int_rules.push_back(ir_itr->second);
360 
361  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases_a = pb_a.getBases();
362  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases_a.begin();
363  b_itr != bases_a.end(); ++b_itr)
364  needs_a.bases.push_back(b_itr->second);
365 
366  WorksetNeeds needs_b;
367  needs_b.cellData = pb_b.cellData();
368 
369  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules_b = pb_b.getIntegrationRules();
370  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules_b.begin();
371  ir_itr != int_rules_b.end(); ++ir_itr)
372  needs_b.int_rules.push_back(ir_itr->second);
373 
374  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases_b = pb_b.getBases();
375  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases_b.begin();
376  b_itr != bases_b.end(); ++b_itr)
377  needs_b.bases.push_back(b_itr->second);
378 
379  return buildEdgeWorksets(needs_a,pb_a.elementBlockID(),local_cell_ids_a,local_side_ids_a,vertex_coordinates_a,
380  needs_b,pb_b.elementBlockID(),local_cell_ids_b,local_side_ids_b,vertex_coordinates_b);
381 /*
382  using std::vector;
383  using std::string;
384  using Teuchos::RCP;
385  using Teuchos::rcp;
386 
387  panzer::MDFieldArrayFactory mdArrayFactory("",true);
388 
389  std::size_t total_num_cells_a = local_cell_ids_a.size();
390  std::size_t total_num_cells_b = local_cell_ids_b.size();
391 
392  TEUCHOS_ASSERT(total_num_cells_a==total_num_cells_b);
393  TEUCHOS_ASSERT(local_side_ids_a.size() == local_cell_ids_a.size());
394  TEUCHOS_ASSERT(local_side_ids_a.size() == static_cast<std::size_t>(vertex_coordinates_a.dimension(0)));
395  TEUCHOS_ASSERT(local_side_ids_b.size() == local_cell_ids_b.size());
396  TEUCHOS_ASSERT(local_side_ids_b.size() == static_cast<std::size_t>(vertex_coordinates_b.dimension(0)));
397 
398  std::size_t total_num_cells = total_num_cells_a;
399 
400  std::size_t workset_size = pb_a.cellData().numCells();
401 
402  Teuchos::RCP< std::vector<panzer::Workset> > worksets_ptr =
403  Teuchos::rcp(new std::vector<panzer::Workset>);
404  std::vector<panzer::Workset>& worksets = *worksets_ptr;
405 
406  // special case for 0 elements!
407  if(total_num_cells==0) {
408 
409  // Setup integration rules and basis
410  RCP<vector<int> > ir_degrees = rcp(new vector<int>(0));
411  RCP<vector<string> > basis_names = rcp(new vector<string>(0));
412 
413  worksets.resize(1);
414  std::vector<panzer::Workset>::iterator i = worksets.begin();
415 
416  i->details(0)->block_id = pb_a.elementBlockID();
417  i->other = Teuchos::rcp(new panzer::WorksetDetails);
418  i->details(1).block_id = pb_b.elementBlockID();
419 
420  i->num_cells = 0;
421  i->ir_degrees = ir_degrees;
422  i->basis_names = basis_names;
423 
424  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules = pb_a.getIntegrationRules();
425 
426  for (std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules.begin();
427  ir_itr != int_rules.end(); ++ir_itr) {
428 
429  RCP<panzer::IntegrationValues2<double> > iv2 =
430  rcp(new panzer::IntegrationValues2<double>("",true));
431  iv2->setupArrays(ir_itr->second);
432 
433  ir_degrees->push_back(ir_itr->first);
434  // i->int_rules.push_back(iv);
435  i->int_rules.push_back(iv2);
436  }
437 
438  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases = pb_a.getBases();
439 
440  // Need to create all combinations of basis/ir pairings
441  for (std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules.begin();
442  ir_itr != int_rules.end(); ++ir_itr) {
443 
444  for (std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases.begin();
445  b_itr != bases.end(); ++b_itr) {
446 
447  RCP<panzer::BasisIRLayout> b_layout = rcp(new panzer::BasisIRLayout(b_itr->second,*ir_itr->second));
448 
449  RCP<panzer::BasisValues2<double> > bv2 =
450  rcp(new panzer::BasisValues2<double>("",true,true));
451  bv2->setupArrays(b_layout);
452  i->bases.push_back(bv2);
453 
454  basis_names->push_back(b_layout->name());
455  }
456 
457  }
458 
459  return worksets_ptr;
460  } // end special case
461 
462  // This collects all the elements that share the same sub cell pairs, this makes it easier to
463  // build the required worksets
464  // key is the pair of local face indices, value is a vector of cell indices that satisfy this pair
465  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > element_list;
466  for (std::size_t cell=0; cell < local_cell_ids_a.size(); ++cell)
467  element_list[std::make_pair(local_side_ids_a[cell],local_side_ids_b[cell])].push_back(cell);
468 
469  // this is the lone iterator that will be used to loop over the element edge list
470  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator edge;
471 
472  // figure out how many worksets will be needed, resize workset vector accordingly
473  std::size_t num_worksets = 0;
474  for(edge=element_list.begin(); edge!=element_list.end();++edge) {
475  std::size_t num_worksets_for_edge = edge->second.size() / workset_size;
476  std::size_t last_workset_size = edge->second.size() % workset_size;
477  if(last_workset_size!=0)
478  num_worksets_for_edge += 1;
479 
480  num_worksets += num_worksets_for_edge;
481  }
482  worksets.resize(num_worksets);
483 
484  // fill the worksets
485  std::vector<Workset>::iterator current_workset = worksets.begin();
486  for(edge=element_list.begin(); edge!=element_list.end();++edge) {
487  // loop over each workset
488  const std::vector<std::size_t> & cell_indices = edge->second;
489 
490  current_workset = buildEdgeWorksets(cell_indices,
491  pb_a,local_cell_ids_a,local_side_ids_a,vertex_coordinates_a,
492  pb_b,local_cell_ids_b,local_side_ids_b,vertex_coordinates_b,
493  current_workset);
494  }
495 
496  // sanity check
497  TEUCHOS_ASSERT(current_workset==worksets.end());
498 
499  return worksets_ptr;
500 */
501 }
502 
503 template<typename ArrayT>
504 Teuchos::RCP<std::vector<panzer::Workset> >
505 panzer::buildEdgeWorksets(const WorksetNeeds & needs_a,
506  const std::string & eblock_a,
507  const std::vector<std::size_t>& local_cell_ids_a,
508  const std::vector<std::size_t>& local_side_ids_a,
509  const ArrayT& vertex_coordinates_a,
510  const WorksetNeeds & needs_b,
511  const std::string & eblock_b,
512  const std::vector<std::size_t>& local_cell_ids_b,
513  const std::vector<std::size_t>& local_side_ids_b,
514  const ArrayT& vertex_coordinates_b)
515 {
516  using Teuchos::RCP;
517  using Teuchos::rcp;
518 
519  panzer::MDFieldArrayFactory mdArrayFactory("",true);
520 
521  std::size_t total_num_cells_a = local_cell_ids_a.size();
522  std::size_t total_num_cells_b = local_cell_ids_b.size();
523 
524  TEUCHOS_ASSERT(total_num_cells_a==total_num_cells_b);
525  TEUCHOS_ASSERT(local_side_ids_a.size() == local_cell_ids_a.size());
526  TEUCHOS_ASSERT(local_side_ids_a.size() == static_cast<std::size_t>(vertex_coordinates_a.dimension(0)));
527  TEUCHOS_ASSERT(local_side_ids_b.size() == local_cell_ids_b.size());
528  TEUCHOS_ASSERT(local_side_ids_b.size() == static_cast<std::size_t>(vertex_coordinates_b.dimension(0)));
529 
530  std::size_t total_num_cells = total_num_cells_a;
531 
532  std::size_t workset_size = needs_a.cellData.numCells();
533 
534  Teuchos::RCP< std::vector<panzer::Workset> > worksets_ptr =
535  Teuchos::rcp(new std::vector<panzer::Workset>);
536  std::vector<panzer::Workset>& worksets = *worksets_ptr;
537 
538  // special case for 0 elements!
539  if(total_num_cells==0) {
540 
541  // Setup integration rules and basis
542  RCP<std::vector<int> > ir_degrees = rcp(new std::vector<int>(0));
543  RCP<std::vector<std::string> > basis_names = rcp(new std::vector<std::string>(0));
544 
545  worksets.resize(1);
546  std::vector<panzer::Workset>::iterator i = worksets.begin();
547 
548  i->details(0).block_id = eblock_a;
549  i->other = Teuchos::rcp(new panzer::WorksetDetails);
550  i->details(1).block_id = eblock_b;
551 
552  i->num_cells = 0;
553  i->ir_degrees = ir_degrees;
554  i->basis_names = basis_names;
555 
556  for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
557 
558  RCP<panzer::IntegrationValues2<double> > iv2 =
559  rcp(new panzer::IntegrationValues2<double>("",true));
560  iv2->setupArrays(needs_a.int_rules[j]);
561 
562  ir_degrees->push_back(needs_a.int_rules[j]->cubature_degree);
563  i->int_rules.push_back(iv2);
564  }
565 
566  // Need to create all combinations of basis/ir pairings
567  for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
568 
569  for(std::size_t b=0;b<needs_a.bases.size();b++) {
570 
571  RCP<panzer::BasisIRLayout> b_layout = rcp(new panzer::BasisIRLayout(needs_a.bases[b],*needs_a.int_rules[j]));
572 
573  RCP<panzer::BasisValues2<double> > bv2 =
574  rcp(new panzer::BasisValues2<double>("",true,true));
575  bv2->setupArrays(b_layout);
576  i->bases.push_back(bv2);
577 
578  basis_names->push_back(b_layout->name());
579  }
580 
581  }
582 
583  return worksets_ptr;
584  } // end special case
585 
586  // This collects all the elements that share the same sub cell pairs, this makes it easier to
587  // build the required worksets
588  // key is the pair of local face indices, value is a vector of cell indices that satisfy this pair
589  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> > element_list;
590  for (std::size_t cell=0; cell < local_cell_ids_a.size(); ++cell)
591  element_list[std::make_pair(local_side_ids_a[cell],local_side_ids_b[cell])].push_back(cell);
592 
593  // this is the lone iterator that will be used to loop over the element edge list
594  std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator edge;
595 
596  // figure out how many worksets will be needed, resize workset vector accordingly
597  std::size_t num_worksets = 0;
598  for(edge=element_list.begin(); edge!=element_list.end();++edge) {
599  std::size_t num_worksets_for_edge = edge->second.size() / workset_size;
600  std::size_t last_workset_size = edge->second.size() % workset_size;
601  if(last_workset_size!=0)
602  num_worksets_for_edge += 1;
603 
604  num_worksets += num_worksets_for_edge;
605  }
606  worksets.resize(num_worksets);
607 
608  // fill the worksets
609  std::vector<Workset>::iterator current_workset = worksets.begin();
610  for(edge=element_list.begin(); edge!=element_list.end();++edge) {
611  // loop over each workset
612  const std::vector<std::size_t> & cell_indices = edge->second;
613 
614  current_workset = buildEdgeWorksets(cell_indices,
615  needs_a,eblock_a,local_cell_ids_a,local_side_ids_a,vertex_coordinates_a,
616  needs_b,eblock_b,local_cell_ids_b,local_side_ids_b,vertex_coordinates_b,
617  current_workset);
618  }
619 
620  // sanity check
621  TEUCHOS_ASSERT(current_workset==worksets.end());
622 
623  return worksets_ptr;
624 }
625 
626 template<typename ArrayT>
627 std::vector<panzer::Workset>::iterator
628 panzer::buildEdgeWorksets(const std::vector<std::size_t> & cell_indices,
629  const WorksetNeeds & needs_a,
630  const std::string & eblock_a,
631  const std::vector<std::size_t>& local_cell_ids_a,
632  const std::vector<std::size_t>& local_side_ids_a,
633  const ArrayT& vertex_coordinates_a,
634  const WorksetNeeds & needs_b,
635  const std::string & eblock_b,
636  const std::vector<std::size_t>& local_cell_ids_b,
637  const std::vector<std::size_t>& local_side_ids_b,
638  const ArrayT& vertex_coordinates_b,
639  std::vector<Workset>::iterator beg)
640 {
641  panzer::MDFieldArrayFactory mdArrayFactory("",true);
642 
643  std::vector<Workset>::iterator wkst = beg;
644 
645  std::size_t current_cell_index = 0;
646  while (current_cell_index<cell_indices.size()) {
647  std::size_t workset_size = needs_a.cellData.numCells();
648 
649  // allocate workset details (details(0) is already associated with the
650  // workset object itself)
651  wkst->other = Teuchos::rcp(new panzer::WorksetDetails);
652 
653  wkst->subcell_dim = needs_a.cellData.baseCellDimension()-1;
654 
655  wkst->details(0).subcell_index = local_side_ids_a[cell_indices[current_cell_index]];
656  wkst->details(0).block_id = eblock_a;
657  wkst->details(0).cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
658  vertex_coordinates_a.dimension(1),
659  vertex_coordinates_a.dimension(2));
660 
661  wkst->details(1).subcell_index = local_side_ids_b[cell_indices[current_cell_index]];
662  wkst->details(1).block_id = eblock_b;
663  wkst->details(1).cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("cvc",workset_size,
664  vertex_coordinates_a.dimension(1),
665  vertex_coordinates_a.dimension(2));
666 
667  std::size_t remaining_cells = cell_indices.size()-current_cell_index;
668  if(remaining_cells<workset_size)
669  workset_size = remaining_cells;
670 
671  // this is the true number of cells in this workset
672  wkst->num_cells = workset_size;
673  wkst->details(0).cell_local_ids.resize(workset_size);
674  wkst->details(1).cell_local_ids.resize(workset_size);
675 
676  for(std::size_t cell=0;cell<workset_size; cell++,current_cell_index++) {
677 
678  wkst->details(0).cell_local_ids[cell] = local_cell_ids_a[cell_indices[current_cell_index]];
679  wkst->details(1).cell_local_ids[cell] = local_cell_ids_b[cell_indices[current_cell_index]];
680 
681  for (std::size_t vertex = 0; vertex < Teuchos::as<std::size_t>(vertex_coordinates_a.dimension(1)); ++ vertex) {
682  for (std::size_t dim = 0; dim < Teuchos::as<std::size_t>(vertex_coordinates_a.dimension(2)); ++ dim) {
683  wkst->details(0).cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates_a(cell_indices[current_cell_index],vertex,dim);
684  wkst->details(1).cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates_b(cell_indices[current_cell_index],vertex,dim);
685  }
686  }
687  }
688 
689  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());
690  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());
691  for(std::size_t i=0;i<wkst->details(0).cell_local_ids.size();i++)
692  cell_local_ids_k_0(i) = wkst->details(0).cell_local_ids[i];
693  for(std::size_t i=0;i<wkst->details(1).cell_local_ids.size();i++)
694  cell_local_ids_k_1(i) = wkst->details(1).cell_local_ids[i];
695  wkst->details(0).cell_local_ids_k = cell_local_ids_k_0;
696  wkst->details(1).cell_local_ids_k = cell_local_ids_k_1;
697 
698  // fill the BasisValues and IntegrationValues arrays
699  std::size_t max_workset_size = needs_a.cellData.numCells();
700  populateValueArrays(max_workset_size,true,needs_a,wkst->details(0)); // populate "side" values
701  populateValueArrays(max_workset_size,true,needs_b,wkst->details(1),Teuchos::rcpFromRef(wkst->details(0)));
702 
703  wkst++;
704  }
705 
706  return wkst;
707 }
708 
709 template<typename ArrayT>
710 std::vector<panzer::Workset>::iterator
711 panzer::buildEdgeWorksets(const std::vector<std::size_t> & cell_indices,
712  const panzer::PhysicsBlock & pb_a,
713  const std::vector<std::size_t>& local_cell_ids_a,
714  const std::vector<std::size_t>& local_side_ids_a,
715  const ArrayT& vertex_coordinates_a,
716  const panzer::PhysicsBlock & pb_b,
717  const std::vector<std::size_t>& local_cell_ids_b,
718  const std::vector<std::size_t>& local_side_ids_b,
719  const ArrayT& vertex_coordinates_b,
720  std::vector<Workset>::iterator beg)
721 {
722  using Teuchos::RCP;
723  using Teuchos::rcp;
724 
725  WorksetNeeds needs_a;
726  needs_a.cellData = pb_a.cellData();
727 
728  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules_a = pb_a.getIntegrationRules();
729  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules_a.begin();
730  ir_itr != int_rules_a.end(); ++ir_itr)
731  needs_a.int_rules.push_back(ir_itr->second);
732 
733  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases_a = pb_a.getBases();
734  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases_a.begin();
735  b_itr != bases_a.end(); ++b_itr)
736  needs_a.bases.push_back(b_itr->second);
737 
738  WorksetNeeds needs_b;
739  needs_b.cellData = pb_b.cellData();
740 
741  const std::map<int,RCP<panzer::IntegrationRule> >& int_rules_b = pb_b.getIntegrationRules();
742  for(std::map<int,RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules_b.begin();
743  ir_itr != int_rules_b.end(); ++ir_itr)
744  needs_b.int_rules.push_back(ir_itr->second);
745 
746  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases_b = pb_b.getBases();
747  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases_b.begin();
748  b_itr != bases_b.end(); ++b_itr)
749  needs_b.bases.push_back(b_itr->second);
750 
751  return buildEdgeWorksets(cell_indices,
752  needs_a,pb_a.elementBlockID(),local_cell_ids_a,local_side_ids_a,vertex_coordinates_a,
753  needs_b,pb_b.elementBlockID(),local_cell_ids_b,local_side_ids_b,vertex_coordinates_b,
754  beg);
755 }
756 
757 namespace panzer {
758 namespace impl {
759 /* Associate two sets of local side IDs into lists. Each list L has the property
760  * that every local side id in that list is the same, and this holds for each
761  * local side ID set. The smallest set of lists is found. The motivation for
762  * this procedure is to find a 1-1 workset pairing in advance. See the comment
763  * re: Intrepid2 in buildBCWorkset for more.
764  * The return value is an RCP to a map. Only the map's values are of interest
765  * in practice. Each value is a list L. The map's key is a pair (side ID a, side
766  * ID b) that gives rise to the list. We return a pointer to a map so that the
767  * caller does not have to deal with the map type; auto suffices.
768  */
769 Teuchos::RCP< std::map<std::pair<std::size_t, std::size_t>, std::vector<std::size_t> > >
770 associateCellsBySideIds(const std::vector<std::size_t>& sia /* local_side_ids_a */,
771  const std::vector<std::size_t>& sib /* local_side_ids_b */)
772 {
773  TEUCHOS_ASSERT(sia.size() == sib.size());
774 
775  auto sip2i_p = Teuchos::rcp(new std::map< std::pair<std::size_t, std::size_t>, std::vector<std::size_t> >);
776  auto& sip2i = *sip2i_p;
777 
778  for (std::size_t i = 0; i < sia.size(); ++i)
779  sip2i[std::make_pair(sia[i], sib[i])].push_back(i);
780 
781  return sip2i_p;
782 }
783 
784 template<typename ArrayT>
785 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
787  const std::vector<std::size_t>& local_cell_ids_a,
788  const std::vector<std::size_t>& local_side_ids_a,
789  const ArrayT& vertex_coordinates_a,
790  const panzer::PhysicsBlock& pb_b,
791  const std::vector<std::size_t>& local_cell_ids_b,
792  const std::vector<std::size_t>& local_side_ids_b,
793  const ArrayT& vertex_coordinates_b,
794  const WorksetNeeds& needs_b)
795 {
796  TEUCHOS_ASSERT(local_cell_ids_a.size() == local_cell_ids_b.size());
797  // Get a and b workset maps separately, but don't populate b's arrays.
798  const Teuchos::RCP<std::map<unsigned,panzer::Workset> >
799  mwa = buildBCWorkset(pb_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a),
800  mwb = buildBCWorkset(needs_b, pb_b.elementBlockID(), local_cell_ids_b, local_side_ids_b,
801  vertex_coordinates_b, false /* populate_value_arrays */);
802  TEUCHOS_ASSERT(mwa->size() == 1 && mwb->size() == 1);
803  for (std::map<unsigned,panzer::Workset>::iterator ait = mwa->begin(), bit = mwb->begin();
804  ait != mwa->end(); ++ait, ++bit) {
805  TEUCHOS_ASSERT(Teuchos::as<std::size_t>(ait->second.num_cells) == local_cell_ids_a.size() &&
806  Teuchos::as<std::size_t>(bit->second.num_cells) == local_cell_ids_b.size());
807  panzer::Workset& wa = ait->second;
808  // Copy b's details(0) to a's details(1).
809  wa.other = Teuchos::rcp(new panzer::WorksetDetails(bit->second.details(0)));
810  // Populate details(1) arrays so that IP are in order corresponding to details(0).
811  populateValueArrays(wa.num_cells, true, needs_b, wa.details(1), Teuchos::rcpFromRef(wa.details(0)));
812  }
813  // Now mwa has everything we need.
814  return mwa;
815 }
816 
817 // Set s = a(idxs). No error checking.
818 template <typename T>
819 void subset(const std::vector<T>& a, const std::vector<std::size_t>& idxs, std::vector<T>& s)
820 {
821  s.resize(idxs.size());
822  for (std::size_t i = 0; i < idxs.size(); ++i)
823  s[i] = a[idxs[i]];
824 }
825 } // namespace impl
826 } // namespace panzer
827 
828 template<typename ArrayT>
829 Teuchos::RCP<std::map<unsigned,panzer::Workset> >
831  const std::vector<std::size_t>& local_cell_ids_a,
832  const std::vector<std::size_t>& local_side_ids_a,
833  const ArrayT& vertex_coordinates_a,
834  const panzer::PhysicsBlock& pb_b,
835  const std::vector<std::size_t>& local_cell_ids_b,
836  const std::vector<std::size_t>& local_side_ids_b,
837  const ArrayT& vertex_coordinates_b)
838 {
839  // Get b's needs.
840  WorksetNeeds needs_b;
841  {
842  needs_b.cellData = pb_b.cellData();
843  const std::map<int,Teuchos::RCP<panzer::IntegrationRule> >& int_rules = pb_b.getIntegrationRules();
844  for (std::map<int,Teuchos::RCP<panzer::IntegrationRule> >::const_iterator ir_itr = int_rules.begin();
845  ir_itr != int_rules.end(); ++ir_itr)
846  needs_b.int_rules.push_back(ir_itr->second);
847  const std::map<std::string,Teuchos::RCP<panzer::PureBasis> >& bases= pb_b.getBases();
848  for(std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator b_itr = bases.begin();
849  b_itr != bases.end(); ++b_itr)
850  needs_b.bases.push_back(b_itr->second);
851  }
852  // Since Intrepid2 requires all side IDs in a workset to be the same (see
853  // Intrepid2 comment above), we break the element list into pieces such that
854  // each piece contains elements on each side of the interface L_a and L_b and
855  // all elemnets L_a have the same side ID, and the same for L_b.
856  auto side_id_associations = impl::associateCellsBySideIds(local_side_ids_a, local_side_ids_b);
857  if (side_id_associations->size() == 1) {
858  // Common case of one workset on each side; optimize for it.
859  return impl::buildBCWorksetForUniqueSideId(pb_a, local_cell_ids_a, local_side_ids_a, vertex_coordinates_a,
860  pb_b, local_cell_ids_b, local_side_ids_b, vertex_coordinates_b,
861  needs_b);
862  } else {
863  // The interface has elements having a mix of side IDs, so deal with each
864  // pair in turn.
865  Teuchos::RCP<std::map<unsigned,panzer::Workset> > mwa = Teuchos::rcp(new std::map<unsigned,panzer::Workset>);
866  std::vector<std::size_t> lci_a, lci_b, lsi_a, lsi_b;
867  panzer::MDFieldArrayFactory mdArrayFactory("", true);
868  const int d1 = Teuchos::as<std::size_t>(vertex_coordinates_a.dimension(1)),
869  d2 = Teuchos::as<std::size_t>(vertex_coordinates_a.dimension(2));
870  for (auto it : *side_id_associations) {
871  const auto& idxs = it.second;
872  impl::subset(local_cell_ids_a, idxs, lci_a);
873  impl::subset(local_side_ids_a, idxs, lsi_a);
874  impl::subset(local_cell_ids_b, idxs, lci_b);
875  impl::subset(local_side_ids_b, idxs, lsi_b);
876  auto vc_a = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("vc_a", idxs.size(), d1, d2);
877  auto vc_b = mdArrayFactory.buildStaticArray<double,Cell,NODE,Dim>("vc_b", idxs.size(), d1, d2);
878  for (std::size_t i = 0; i < idxs.size(); ++i) {
879  const auto ii = idxs[i];
880  for (int j = 0; j < d1; ++j)
881  for (int k = 0; k < d2; ++k) {
882  vc_a(i, j, k) = vertex_coordinates_a(ii, j, k);
883  vc_b(i, j, k) = vertex_coordinates_b(ii, j, k);
884  }
885  }
886  auto mwa_it = impl::buildBCWorksetForUniqueSideId(pb_a, lci_a, lsi_a, vc_a, pb_b, lci_b, lsi_b, vc_b, needs_b);
887  TEUCHOS_ASSERT(mwa_it->size() == 1);
888  // Form a unique key that encodes the pair (side ID a, side ID b). We
889  // abuse the key here in the sense that it is everywhere else understood
890  // to correspond to the side ID of the elements in the workset.
891  // 1000 is a number substantially larger than is needed for any element.
892  const std::size_t key = lsi_a[0] * 1000 + lsi_b[0];
893  (*mwa)[key] = mwa_it->begin()->second;
894  }
895  return mwa;
896  }
897 }
898 
899 #endif
Teuchos::RCP< std::map< unsigned, Workset > > buildBCWorkset(const PhysicsBlock &volume_pb, const std::vector< std::size_t > &local_cell_ids, const std::vector< std::size_t > &local_side_ids, const ArrayT &vertex_coordinates)
Object that contains information on the physics and discretization of a block of elements with the SA...
std::string elementBlockID() const
PHX::MDField< ScalarT > vector
Teuchos::RCP< std::vector< Workset > > buildEdgeWorksets(const PhysicsBlock &pb_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 PhysicsBlock &pb_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::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< std::vector< Workset > > buildWorksets(const PhysicsBlock &pb, const std::vector< std::size_t > &local_cell_ids, const ArrayT &vertex_coordinates)
const panzer::CellData & cellData() const
void subset(const std::vector< T > &a, const std::vector< std::size_t > &idxs, std::vector< T > &s)
const std::map< std::string, Teuchos::RCP< panzer::PureBasis > > & getBases() const
Returns the unique set of bases, key is the unique panzer::PureBasis::name() of the basis...
Teuchos::RCP< WorksetDetails > other
void populateValueArrays(std::size_t num_cells, bool isSide, const WorksetNeeds &needs, WorksetDetails &details, const Teuchos::RCP< WorksetDetails > other_details)
const std::map< int, Teuchos::RCP< panzer::IntegrationRule > > & getIntegrationRules() const
Returns the unique set of point rules, key is the unique panzer::PointRule::name() ...
Teuchos::RCP< std::map< unsigned, panzer::Workset > > buildBCWorksetForUniqueSideId(const panzer::PhysicsBlock &pb_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::PhysicsBlock &pb_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_b)