43 #ifndef __Panzer_Workset_Builder_impl_hpp__
44 #define __Panzer_Workset_Builder_impl_hpp__
56 #include "Phalanx_DataLayout_MDALayout.hpp"
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"
65 template<
typename ArrayT>
68 const std::string & elementBlock,
69 const std::vector<std::size_t>& local_cell_ids,
70 const ArrayT& vertex_coordinates)
79 std::size_t total_num_cells = local_cell_ids.size();
81 std::size_t workset_size = needs.cellData.numCells();
85 std::vector<panzer::Workset>& worksets = *worksets_ptr;
88 if(total_num_cells==0) {
91 RCP<vector<int> > ir_degrees =
rcp(
new vector<int>(0));
92 RCP<vector<string> > basis_names =
rcp(
new vector<string>(0));
95 std::vector<panzer::Workset>::iterator i = worksets.begin();
97 i->block_id = elementBlock;
98 i->ir_degrees = ir_degrees;
99 i->basis_names = basis_names;
101 for (std::size_t j=0;j<needs.int_rules.size();j++) {
103 RCP<panzer::IntegrationValues2<double> > iv2 =
105 iv2->setupArrays(needs.int_rules[j]);
107 ir_degrees->push_back(needs.int_rules[j]->cubature_degree);
108 i->int_rules.push_back(iv2);
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
117 RCP<panzer::BasisValues2<double> > bv2
119 bv2->setupArrays(b_layout);
120 i->bases.push_back(bv2);
122 basis_names->push_back(b_layout->name());
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) {
136 last_set_is_full =
false;
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;
144 if (!last_set_is_full) {
145 worksets.back().num_cells = last_workset_size;
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);
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;
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;
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) {
179 wkst->cell_vertex_coordinates(cell,vertex,dim) = vertex_coordinates(cell + offset,vertex,dim);
182 offset += wkst->num_cells;
185 TEUCHOS_ASSERT(offset == Teuchos::as<std::size_t>(vertex_coordinates.extent(0)));
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;
196 for(std::vector<panzer::Workset>::iterator wkst = worksets.begin(); wkst != worksets.end(); ++wkst)
205 template<
typename ArrayT>
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)
231 TEUCHOS_ASSERT(local_side_ids.size() ==
static_cast<std::size_t
>(vertex_coordinates.extent(0)));
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]));
238 std::map<unsigned,panzer::Workset>& worksets = *worksets_ptr;
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) {
244 std::vector<std::size_t>& cell_local_ids = worksets[side->first].cell_local_ids;
246 worksets[side->first].cell_vertex_coordinates = mdArrayFactory.buildStaticArray<double,Cell,
NODE,Dim>(
"cvc",
248 vertex_coordinates.extent(1),
249 vertex_coordinates.extent(2));
250 Workset::CellCoordArray coords = worksets[side->first].cell_vertex_coordinates;
252 for (std::size_t cell = 0; cell < side->second.size(); ++cell) {
253 cell_local_ids.push_back(side->second[cell].second);
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);
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;
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;
272 if (populate_value_arrays) {
274 for (std::map<unsigned,panzer::Workset>::iterator wkst = worksets.begin();
275 wkst != worksets.end(); ++wkst) {
302 const std::vector<std::size_t>& sib )
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;
309 for (std::size_t i = 0; i < sia.size(); ++i)
310 sip2i[std::make_pair(sia[i], sib[i])].push_back(i);
316 template <
typename T>
317 void subset(
const std::vector<T>& a,
const std::vector<std::size_t>& idxs, std::vector<T>& s)
319 s.resize(idxs.size());
320 for (std::size_t i = 0; i < idxs.size(); ++i)
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,
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,
338 TEUCHOS_ASSERT(local_cell_ids_a.size() == local_cell_ids_b.size());
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 );
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());
353 populateValueArrays(wa.num_cells,
true, needs_b, wa.details(1), Teuchos::rcpFromRef(wa.details(0)));
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,
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)
383 if (side_id_associations->size() == 1) {
386 needs_b, blockid_b, local_cell_ids_b, local_side_ids_b, vertex_coordinates_b,
392 std::vector<std::size_t> lci_a, lci_b, lsi_a, lsi_b;
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;
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);
413 needs_b,blockid_b, lci_b, lsi_b, vc_b,
420 const std::size_t key = lsi_a[0] * 1000 + lsi_b[0];
421 (*mwa)[key] = mwa_it->begin()->second;
430 template<
typename ArrayT>
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)
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();
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)));
457 std::size_t total_num_cells = total_num_cells_a;
459 std::size_t workset_size = needs_a.cellData.numCells();
463 std::vector<panzer::Workset>& worksets = *worksets_ptr;
466 if(total_num_cells==0) {
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));
473 std::vector<panzer::Workset>::iterator i = worksets.begin();
475 i->details(0).block_id = eblock_a;
477 i->details(1).block_id = eblock_b;
480 i->ir_degrees = ir_degrees;
481 i->basis_names = basis_names;
483 for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
485 RCP<panzer::IntegrationValues2<double> > iv2 =
487 iv2->setupArrays(needs_a.int_rules[j]);
489 ir_degrees->push_back(needs_a.int_rules[j]->cubature_degree);
490 i->int_rules.push_back(iv2);
494 for(std::size_t j=0;j<needs_a.int_rules.size();j++) {
496 for(std::size_t b=0;b<needs_a.bases.size();b++) {
500 RCP<panzer::BasisValues2<double> > bv2 =
502 bv2->setupArrays(b_layout);
503 i->bases.push_back(bv2);
505 basis_names->push_back(b_layout->name());
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);
521 std::map<std::pair<unsigned,unsigned>,std::vector<std::size_t> >::const_iterator edge;
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;
531 num_worksets += num_worksets_for_edge;
533 worksets.resize(num_worksets);
536 std::vector<Workset>::iterator current_workset = worksets.begin();
537 for(edge=element_list.begin(); edge!=element_list.end();++edge) {
539 const std::vector<std::size_t> & cell_indices = edge->second;
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,
553 template<
typename ArrayT>
554 std::vector<panzer::Workset>::iterator
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)
570 std::vector<Workset>::iterator wkst = beg;
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();
580 wkst->subcell_dim = needs_a.cellData.baseCellDimension()-1;
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));
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));
594 std::size_t remaining_cells = cell_indices.size()-current_cell_index;
595 if(remaining_cells<workset_size)
596 workset_size = remaining_cells;
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);
603 for(std::size_t cell=0;cell<workset_size; cell++,current_cell_index++) {
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]];
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);
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;
626 std::size_t max_workset_size = needs_a.cellData.numCells();
628 populateValueArrays(max_workset_size,
true,needs_b,wkst->details(1),Teuchos::rcpFromRef(wkst->details(0)));
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)