Panzer Version of the Day
Loading...
Searching...
No Matches
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
65template<typename ArrayT>
66Teuchos::RCP< std::vector<panzer::Workset> >
67panzer::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
209template<typename ArrayT>
210Teuchos::RCP<std::map<unsigned,panzer::Workset> >
211panzer::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
305namespace panzer {
306namespace 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 */
318Teuchos::RCP< std::map<std::pair<std::size_t, std::size_t>, std::vector<std::size_t> > >
319associateCellsBySideIds(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.
334template <typename T>
335void 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
342template<typename ArrayT>
343Teuchos::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
383template<typename ArrayT>
384Teuchos::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
456template<typename ArrayT>
457Teuchos::RCP<std::vector<panzer::Workset> >
458panzer::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
576template<typename ArrayT>
577std::vector<panzer::Workset>::iterator
578panzer::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
int num_cells
DEPRECATED - use: numCells()
WorksetDetails & details(const int i)
Convenience wrapper to operator() for pointer access.
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)
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)
void populateValueArrays(std::size_t num_cells, bool isSide, const WorksetNeeds &needs, WorksetDetails &details, const Teuchos::RCP< WorksetDetails > other_details)
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::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)