Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_L2Projection_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // @HEADER
3 
4 #ifndef PANZER_L2_PROJECTION_IMPL_HPP
5 #define PANZER_L2_PROJECTION_IMPL_HPP
6 
7 #include "Teuchos_DefaultMpiComm.hpp"
8 #include "Tpetra_CrsGraph.hpp"
9 #include "Tpetra_MultiVector.hpp"
10 #include "Tpetra_CrsMatrix.hpp"
11 #include "Panzer_BasisDescriptor.hpp"
12 #include "Panzer_TpetraLinearObjContainer.hpp"
13 #include "Panzer_TpetraLinearObjFactory.hpp"
14 #include "Panzer_BlockedTpetraLinearObjFactory.hpp"
15 #include "Panzer_ElementBlockIdToPhysicsIdMap.hpp"
16 #include "Panzer_DOFManagerFactory.hpp"
17 #include "Panzer_BlockedDOFManagerFactory.hpp"
18 #include "Panzer_IntrepidFieldPattern.hpp"
19 #include "Panzer_WorksetDescriptor.hpp"
20 #include "Panzer_WorksetContainer.hpp"
21 #include "Panzer_Workset.hpp"
22 
23 namespace panzer {
24 
25  template<typename LO, typename GO>
26  void panzer::L2Projection<LO,GO>::setup(const panzer::BasisDescriptor& targetBasis,
27  const panzer::IntegrationDescriptor& integrationDescriptor,
28  const Teuchos::RCP<const Teuchos::MpiComm<int>>& comm,
29  const Teuchos::RCP<const panzer::ConnManager<LO,GO>>& connManager,
30  const std::vector<std::string>& elementBlockNames,
31  const Teuchos::RCP<panzer::WorksetContainer> worksetContainer)
32  {
33  targetBasisDescriptor_ = targetBasis;
34  integrationDescriptor_ = integrationDescriptor;
35  comm_ = comm;
36  connManager_ = connManager;
37  elementBlockNames_ = elementBlockNames;
38  worksetContainer_ = worksetContainer;
39  setupCalled_ = true;
40 
41  // Build target DOF Manager
42  targetGlobalIndexer_ =
43  Teuchos::rcp(new panzer::DOFManager<LO,GO>(Teuchos::rcp_const_cast<panzer::ConnManager<LO,GO>>(connManager),*(comm->getRawMpiComm())));
44 
45  // For hybrid mesh, blocks could have different topology
46  for (const auto& eBlock : elementBlockNames_) {
47  std::vector<shards::CellTopology> topologies;
48  connManager_->getElementBlockTopologies(topologies);
49  std::vector<std::string> ebNames;
50  connManager_->getElementBlockIds(ebNames);
51  const auto search = std::find(ebNames.cbegin(),ebNames.cend(),eBlock);
52  TEUCHOS_ASSERT(search != ebNames.cend());
53  const int index = std::distance(ebNames.cbegin(),search);
54  const auto& cellTopology = topologies[index];
55 
56  auto intrepidBasis = panzer::createIntrepid2Basis<PHX::Device,double,double>(targetBasisDescriptor_.getType(),
57  targetBasisDescriptor_.getOrder(),
58  cellTopology);
59  Teuchos::RCP<const panzer::FieldPattern> fieldPattern(new panzer::Intrepid2FieldPattern(intrepidBasis));
60  // Field name is the basis type
61  targetGlobalIndexer_->addField(eBlock,targetBasisDescriptor_.getType(),fieldPattern);
62  }
63 
64  targetGlobalIndexer_->buildGlobalUnknowns();
65 
66  // Check workset needs are correct
67  }
68 
69  template<typename LO, typename GO>
70  Teuchos::RCP<panzer::UniqueGlobalIndexer<LO,GO>>
71  panzer::L2Projection<LO,GO>::getTargetGlobalIndexer() const
72  {return targetGlobalIndexer_;}
73 
74  template<typename LO, typename GO>
75  Teuchos::RCP<Tpetra::CrsMatrix<double,LO,GO,panzer::TpetraNodeType>>
76  panzer::L2Projection<LO,GO>::buildMassMatrix(bool use_lumping,
77  const std::unordered_map<std::string,double>* elementBlockMultipliers)
78  {
79  PANZER_FUNC_TIME_MONITOR_DIFF("L2Projection::Build Mass Matrix",TopBuildMassMatrix);
80  TEUCHOS_ASSERT(setupCalled_);
81 
82  if (elementBlockMultipliers != nullptr) {
83  TEUCHOS_ASSERT(elementBlockMultipliers->size() == elementBlockNames_.size());
84  }
85 
86  // Allocate the owned matrix
87  std::vector<Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO>>> indexers;
88  indexers.push_back(targetGlobalIndexer_);
89 
90  panzer::BlockedTpetraLinearObjFactory<panzer::Traits,double,LO,GO,panzer::TpetraNodeType> factory(comm_,indexers);
91 
92  auto ownedMatrix = factory.getTpetraMatrix(0,0);
93  auto ghostedMatrix = factory.getGhostedTpetraMatrix(0,0);
94  ownedMatrix->resumeFill();
95  ownedMatrix->setAllToScalar(0.0);
96  ghostedMatrix->resumeFill();
97  ghostedMatrix->setAllToScalar(0.0);
98  PHX::Device::fence();
99 
100  auto M = ghostedMatrix->getLocalMatrix();
101  const int fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
102 
103  const bool is_scalar = targetBasisDescriptor_.getType()=="HGrad" || targetBasisDescriptor_.getType()=="Const" || targetBasisDescriptor_.getType()=="HVol";
104 
105  // Loop over element blocks and fill mass matrix
106  if(is_scalar){
107  for (const auto& block : elementBlockNames_) {
108 
109  double ebMultiplier = 1.0;
110  if (elementBlockMultipliers != nullptr)
111  ebMultiplier = elementBlockMultipliers->find(block)->second;
112 
113  // Based on descriptor, currently assumes there should only be one workset
114  panzer::WorksetDescriptor wd(block,panzer::WorksetSizeType::ALL_ELEMENTS,true,true);
115  const auto& worksets = worksetContainer_->getWorksets(wd);
116 
117  for (const auto& workset : *worksets) {
118 
119  const auto& basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
120 
121  const auto& unweightedBasis = basisValues.basis_scalar;
122  const auto& weightedBasis = basisValues.weighted_basis_scalar;
123 
124  // Offsets (this assumes UVM, need to fix)
125  const std::vector<LO>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
126  PHX::View<LO*> kOffsets("MassMatrix: Offsets",offsets.size());
127  for(const auto& i : offsets)
128  kOffsets(i) = offsets[i];
129 
130  PHX::Device::fence();
131 
132  // Local Ids
133  Kokkos::View<LO**,PHX::Device> localIds("MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
134  targetGlobalIndexer_->getElementBlockGIDCount(block));
135 
136  // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
137  const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
138 
139  targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
140 
141  const int numBasisPoints = static_cast<int>(weightedBasis.extent(1));
142  if ( use_lumping ) {
143  Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (const int& cell) {
144  double total_mass = 0.0, trace = 0.0;
145 
146  LO cLIDs[256];
147  const int numIds = static_cast<int>(localIds.extent(1));
148  for(int i=0;i<numIds;++i)
149  cLIDs[i] = localIds(cell,i);
150 
151  double vals[256]={0.0};
152  const int numQP = static_cast<int>(unweightedBasis.extent(2));
153 
154  for (int row=0; row < numBasisPoints; ++row) {
155  for (int col=0; col < numIds; ++col) {
156  for (int qp=0; qp < numQP; ++qp) {
157  auto tmp = unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
158  total_mass += tmp;
159  if (col == row )
160  trace += tmp;
161  }
162  }
163  }
164 
165  for (int row=0; row < numBasisPoints; ++row) {
166  for (int col=0; col < numBasisPoints; ++col)
167  vals[col] = 0;
168 
169  int offset = kOffsets(row);
170  LO lid = localIds(cell,offset);
171  int col = row;
172  vals[col] = 0.0;
173  for (int qp=0; qp < numQP; ++qp)
174  vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier * total_mass / trace;
175 
176  M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
177  }
178  });
179 
180  } else {
181  Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (const int& cell) {
182 
183  LO cLIDs[256];
184  const int numIds = static_cast<int>(localIds.extent(1));
185  for(int i=0;i<numIds;++i)
186  cLIDs[i] = localIds(cell,i);
187 
188  double vals[256];
189  const int numQP = static_cast<int>(unweightedBasis.extent(2));
190 
191  for (int row=0; row < numBasisPoints; ++row) {
192  int offset = kOffsets(row);
193  LO lid = localIds(cell,offset);
194 
195  for (int col=0; col < numIds; ++col) {
196  vals[col] = 0.0;
197  for (int qp=0; qp < numQP; ++qp)
198  vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
199  }
200  M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
201 
202  }
203 
204  });
205  }
206  }
207  }
208  } else {
209  for (const auto& block : elementBlockNames_) {
210 
211  double ebMultiplier = 1.0;
212  if (elementBlockMultipliers != nullptr)
213  ebMultiplier = elementBlockMultipliers->find(block)->second;
214 
215  // Based on descriptor, currently assumes there should only be one workset
216  panzer::WorksetDescriptor wd(block,panzer::WorksetSizeType::ALL_ELEMENTS,true,true);
217  const auto& worksets = worksetContainer_->getWorksets(wd);
218 
219  for (const auto& workset : *worksets) {
220 
221  const auto& basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
222 
223  const auto& unweightedBasis = basisValues.basis_vector;
224  const auto& weightedBasis = basisValues.weighted_basis_vector;
225 
226  // Offsets (this assumes UVM, need to fix)
227  const std::vector<LO>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
228  PHX::View<LO*> kOffsets("MassMatrix: Offsets",offsets.size());
229  for(const auto& i : offsets)
230  kOffsets(i) = offsets[i];
231 
232  PHX::Device::fence();
233 
234  // Local Ids
235  Kokkos::View<LO**,PHX::Device> localIds("MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
236  targetGlobalIndexer_->getElementBlockGIDCount(block));
237 
238  // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
239  const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
240 
241  targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
242 
243  const int numBasisPoints = static_cast<int>(weightedBasis.extent(1));
244  Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (const int& cell) {
245 
246  LO cLIDs[256];
247  const int numIds = static_cast<int>(localIds.extent(1));
248  for(int i=0;i<numIds;++i)
249  cLIDs[i] = localIds(cell,i);
250 
251  double vals[256];
252  const int numQP = static_cast<int>(unweightedBasis.extent(2));
253 
254  for (int qp=0; qp < numQP; ++qp) {
255  for (int row=0; row < numBasisPoints; ++row) {
256  int offset = kOffsets(row);
257  LO lid = localIds(cell,offset);
258 
259  for (int col=0; col < numIds; ++col){
260  vals[col] = 0.0;
261  for(int dim=0; dim < static_cast<int>(weightedBasis.extent(3)); ++dim)
262  vals[col] += unweightedBasis(cell,row,qp,dim) * weightedBasis(cell,col,qp,dim) * ebMultiplier;
263  }
264 
265  M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
266  }
267  }
268 
269  });
270  }
271  }
272  }
273  PHX::exec_space::fence();
274 
275  {
276  PANZER_FUNC_TIME_MONITOR_DIFF("Exporting of mass matrix",ExportMM);
277  auto map = factory.getMap(0);
278  ghostedMatrix->fillComplete(map,map);
279  const auto exporter = factory.getGhostedExport(0);
280  ownedMatrix->doExport(*ghostedMatrix, *exporter, Tpetra::ADD);
281  ownedMatrix->fillComplete();
282  }
283  return ownedMatrix;
284  }
285 
286  template<typename LO, typename GO>
287  Teuchos::RCP<Tpetra::MultiVector<double,LO,GO,panzer::TpetraNodeType>>
288  panzer::L2Projection<LO,GO>::buildInverseLumpedMassMatrix()
289  {
290  PANZER_FUNC_TIME_MONITOR_DIFF("L2Projection<LO,GO>::buildInverseLumpedMassMatrix",buildInvLMM);
291  using Teuchos::rcp;
292  const auto massMatrix = this->buildMassMatrix(true);
293  const auto lumpedMassMatrix = rcp(new Tpetra::MultiVector<double,LO,GO,panzer::TpetraNodeType>(massMatrix->getDomainMap(),1,true));
294  const auto tmp = rcp(new Tpetra::MultiVector<double,LO,GO,panzer::TpetraNodeType>(massMatrix->getRangeMap(),1,false));
295  tmp->putScalar(1.0);
296  {
297  PANZER_FUNC_TIME_MONITOR_DIFF("Apply",Apply);
298  massMatrix->apply(*tmp,*lumpedMassMatrix);
299  }
300  {
301  PANZER_FUNC_TIME_MONITOR_DIFF("Reciprocal",Reciprocal);
302  lumpedMassMatrix->reciprocal(*lumpedMassMatrix);
303  }
304  return lumpedMassMatrix;
305  }
306 
307  template<typename LO, typename GO>
308  Teuchos::RCP<Tpetra::CrsMatrix<double,LO,GO,panzer::TpetraNodeType>>
309  panzer::L2Projection<LO,GO>::buildRHSMatrix(const panzer::UniqueGlobalIndexer<LO,GO>& sourceGlobalIndexer,
310  const Teuchos::RCP<const Tpetra::Map<LO,GO,panzer::TpetraNodeType>>& inputOwnedSourceMap,
311  const std::string& sourceFieldName,
312  const panzer::BasisDescriptor& sourceBasisDescriptor,
313  const int directionIndex)
314  {
315  // *******************
316  // Create Retangular matrix (both ghosted and owned).
317  // *******************
318  using Teuchos::RCP;
319  using Teuchos::rcp;
320  using MapType = Tpetra::Map<LO,GO,panzer::TpetraNodeType>;
321  using GraphType = Tpetra::CrsGraph<LO,GO,panzer::TpetraNodeType>;
322  using ExportType = Tpetra::Export<LO,GO,panzer::TpetraNodeType>;
323  using MatrixType = Tpetra::CrsMatrix<double,LO,GO,panzer::TpetraNodeType>;
324 
325  // *******************
326  // Build ghosted graph
327  // *******************
328 
329  RCP<MapType> ghostedTargetMap;
330  {
331  std::vector<GO> indices;
332  targetGlobalIndexer_->getOwnedAndGhostedIndices(indices);
333  ghostedTargetMap = rcp(new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
334  }
335 
336  RCP<MapType> ghostedSourceMap;
337  {
338  std::vector<GO> indices;
339  sourceGlobalIndexer.getOwnedAndGhostedIndices(indices);
340  ghostedSourceMap = rcp(new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
341  }
342 
343  RCP<GraphType> ghostedGraph = rcp(new GraphType(ghostedTargetMap,ghostedSourceMap,0));
344 
345  // Now insert the non-zero pattern per row
346  std::vector<std::string> elementBlockIds;
347  targetGlobalIndexer_->getElementBlockIds(elementBlockIds);
348  std::vector<std::string>::const_iterator blockItr;
349  for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
350  std::string blockId = *blockItr;
351  const std::vector<LO> & elements = targetGlobalIndexer_->getElementBlock(blockId);
352 
353  std::vector<GO> row_gids;
354  std::vector<GO> col_gids;
355 
356  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
357  targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
358  sourceGlobalIndexer.getElementGIDs(elements[elmt],col_gids);
359  for(std::size_t row=0;row<row_gids.size();row++)
360  ghostedGraph->insertGlobalIndices(row_gids[row],col_gids);
361  }
362  }
363 
364  RCP<MapType> ownedTargetMap;
365  {
366  std::vector<GO> indices;
367  targetGlobalIndexer_->getOwnedIndices(indices);
368  ownedTargetMap = rcp(new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
369  }
370 
371  RCP<const MapType> ownedSourceMap = inputOwnedSourceMap;
372  if (is_null(ownedSourceMap)) {
373  std::vector<GO> indices;
374  sourceGlobalIndexer.getOwnedIndices(indices);
375  ownedSourceMap = rcp(new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
376  }
377 
378  // Fill complete with owned range and domain map
379  ghostedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
380 
381  // *****************
382  // Build owned graph
383  // *****************
384 
385  RCP<GraphType> ownedGraph = rcp(new GraphType(ownedTargetMap,0));
386  RCP<const ExportType> exporter = rcp(new ExportType(ghostedTargetMap,ownedTargetMap));
387  ownedGraph->doExport(*ghostedGraph, *exporter, Tpetra::INSERT);
388  ownedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
389 
390  RCP<MatrixType> ghostedMatrix = rcp(new MatrixType(ghostedGraph));
391  RCP<MatrixType> ownedMatrix = rcp(new MatrixType(ownedGraph));
392  // ghostedMatrix.fillComplete();
393  // ghostedMatrix.resumeFill();
394 
395  ghostedMatrix->setAllToScalar(0.0);
396  ownedMatrix->setAllToScalar(0.0);
397  PHX::Device::fence();
398 
399  // *******************
400  // Fill ghosted matrix
401  // *******************
402  for (const auto& block : elementBlockNames_) {
403 
404  panzer::WorksetDescriptor wd(block,panzer::WorksetSizeType::ALL_ELEMENTS,true,true);
405  const auto& worksets = worksetContainer_->getWorksets(wd);
406  for (const auto& workset : *worksets) {
407 
408  // Get target basis values: current implementation assumes target basis is HGrad
409  const auto& targetBasisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
410  const auto& targetWeightedBasis = targetBasisValues.weighted_basis_scalar.get_static_view();
411 
412  // Sources can be any basis
413  const auto& sourceBasisValues = workset.getBasisValues(sourceBasisDescriptor,integrationDescriptor_);
414  Kokkos::View<const double***,PHX::Device> sourceUnweightedScalarBasis;
415  Kokkos::View<const double****,PHX::Device> sourceUnweightedVectorBasis;
416  bool useRankThreeBasis = false; // default to gradient or vector basis
417  if ( (sourceBasisDescriptor.getType() == "HGrad") || (sourceBasisDescriptor.getType() == "Const") || (sourceBasisDescriptor.getType() == "HVol") ) {
418  if (directionIndex == -1) { // Project dof value
419  sourceUnweightedScalarBasis = sourceBasisValues.basis_scalar.get_static_view();
420  useRankThreeBasis = true;
421  }
422  else { // Project dof gradient of scalar basis
423  sourceUnweightedVectorBasis = sourceBasisValues.grad_basis.get_static_view();
424  }
425  }
426  else { // Project vector value
427  sourceUnweightedVectorBasis = sourceBasisValues.basis_vector.get_static_view();
428  }
429 
430  // Get the element local ids
431  Kokkos::View<LO**,PHX::Device> targetLocalIds("buildRHSMatrix: targetLocalIds", workset.numOwnedCells(),
432  targetGlobalIndexer_->getElementBlockGIDCount(block));
433  Kokkos::View<LO**,PHX::Device> sourceLocalIds("buildRHSMatrix: sourceLocalIds", workset.numOwnedCells(),
434  sourceGlobalIndexer.getElementBlockGIDCount(block));
435  {
436  // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
437  const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
438  targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,targetLocalIds);
439  sourceGlobalIndexer.getElementLIDs(cellLocalIdsNoGhost,sourceLocalIds);
440  }
441 
442  // Get the offsets
443  Kokkos::View<LO*,PHX::Device> targetFieldOffsets;
444  {
445  const auto fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
446  const std::vector<LO>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
447  targetFieldOffsets = Kokkos::View<LO*,PHX::Device>("L2Projection:buildRHS:targetFieldOffsets",offsets.size());
448  const auto hostOffsets = Kokkos::create_mirror_view(targetFieldOffsets);
449  for(size_t i=0; i < offsets.size(); ++i)
450  hostOffsets(i) = offsets[i];
451  Kokkos::deep_copy(targetFieldOffsets,hostOffsets);
452  PHX::Device::fence();
453  }
454 
455  Kokkos::View<LO*,PHX::Device> sourceFieldOffsets;
456  {
457  const auto fieldIndex = sourceGlobalIndexer.getFieldNum(sourceFieldName);
458  const std::vector<LO>& offsets = sourceGlobalIndexer.getGIDFieldOffsets(block,fieldIndex);
459  TEUCHOS_TEST_FOR_EXCEPTION(offsets.size() == 0, std::runtime_error,
460  "ERROR: panzer::L2Projection::buildRHSMatrix() - The source field, \""
461  << sourceFieldName << "\", does not exist in element block \""
462  << block << "\"!");
463  sourceFieldOffsets = Kokkos::View<LO*,PHX::Device>("L2Projection:buildRHS:sourceFieldOffsets",offsets.size());
464  const auto hostOffsets = Kokkos::create_mirror_view(sourceFieldOffsets);
465  for(size_t i=0; i <offsets.size(); ++i)
466  hostOffsets(i) = offsets[i];
467  Kokkos::deep_copy(sourceFieldOffsets,hostOffsets);
468  PHX::Device::fence();
469  }
470 
471  const auto localMatrix = ghostedMatrix->getLocalMatrix();
472  const int numRows = static_cast<int>(targetWeightedBasis.extent(1));
473  int tmpNumCols = -1;
474  int tmpNumQP = -1;
475  if (useRankThreeBasis) {
476  tmpNumCols = static_cast<int>(sourceUnweightedScalarBasis.extent(1));
477  tmpNumQP = static_cast<int>(sourceUnweightedScalarBasis.extent(2));
478  }
479  else {
480  tmpNumCols = static_cast<int>(sourceUnweightedVectorBasis.extent(1));
481  tmpNumQP = static_cast<int>(sourceUnweightedVectorBasis.extent(2));
482  }
483  const int numCols = tmpNumCols;
484  const int numQP = tmpNumQP;
485 
486  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.numOwnedCells()),KOKKOS_LAMBDA (const int& cell) {
487  LO cLIDs[256];
488  double vals[256];
489  for (int row = 0; row < numRows; ++row) {
490  const int rowOffset = targetFieldOffsets(row);
491  const int rowLID = targetLocalIds(cell,rowOffset);
492  for (int col = 0; col < numCols; ++col)
493  vals[col] = 0.0;
494 
495  for (int col = 0; col < numCols; ++col) {
496  for (int qp = 0; qp < numQP; ++qp) {
497  const int colOffset = sourceFieldOffsets(col);
498  const int colLID = sourceLocalIds(cell,colOffset);
499  cLIDs[col] = colLID;
500  if (useRankThreeBasis)
501  vals[col] += sourceUnweightedScalarBasis(cell,col,qp) * targetWeightedBasis(cell,row,qp);
502  else
503  vals[col] += sourceUnweightedVectorBasis(cell,col,qp,directionIndex) * targetWeightedBasis(cell,row,qp);
504  }
505  }
506  localMatrix.sumIntoValues(rowLID,cLIDs,numCols,vals,true,true);
507  }
508  });
509  }
510 
511  }
512  ghostedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
513  ownedMatrix->resumeFill();
514  ownedMatrix->doExport(*ghostedMatrix,*exporter,Tpetra::ADD);
515  ownedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
516 
517  return ownedMatrix;
518  }
519 
520 }
521 
522 #endif
panzer::UniqueGlobalIndexer::getElementLIDs
const Kokkos::View< const LocalOrdinalT *, Kokkos::LayoutRight, PHX::Device > getElementLIDs(LocalOrdinalT localElmtId) const
Definition: Panzer_UniqueGlobalIndexer.hpp:302
panzer::UniqueGlobalIndexer::getElementGIDs
virtual void getElementGIDs(LocalOrdinalT localElmtId, std::vector< GlobalOrdinalT > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
panzer::BlockedTpetraLinearObjFactory::getTpetraMatrix
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
Definition: Panzer_BlockedTpetraLinearObjFactory_impl.hpp:1049
Panzer_Workset.hpp
is_null
bool is_null(const boost::shared_ptr< T > &p)
panzer::BlockedTpetraLinearObjFactory::getGhostedExport
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a &quot;normal&quot; object
Definition: Panzer_BlockedTpetraLinearObjFactory_impl.hpp:931
panzer::UniqueGlobalIndexer::getOwnedAndGhostedIndices
virtual void getOwnedAndGhostedIndices(std::vector< GlobalOrdinalT > &indices) const =0
Get the set of owned and ghosted indices for this processor.
panzer::L2Projection::buildInverseLumpedMassMatrix
Teuchos::RCP< Tpetra::MultiVector< double, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildInverseLumpedMassMatrix()
Allocates, fills and returns a Tpetra::MultiVector containing the inverse lumped mass matrix values a...
Definition: Panzer_L2Projection_impl.hpp:288
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
panzer::UniqueGlobalIndexer::getGIDFieldOffsets
virtual const std::vector< int > & getGIDFieldOffsets(const std::string &blockId, int fieldNum) const =0
Use the field pattern so that you can find a particular field in the GIDs array.
Panzer_TpetraLinearObjContainer.hpp
panzer::L2Projection::getTargetGlobalIndexer
Teuchos::RCP< panzer::UniqueGlobalIndexer< LO, GO > > getTargetGlobalIndexer() const
Returns the target global indexer. Will be null if setup() has not been called.
Definition: Panzer_L2Projection_impl.hpp:71
Panzer_BasisDescriptor.hpp
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
panzer::BasisDescriptor::getType
const std::string & getType() const
Get type of basis.
Definition: Panzer_BasisDescriptor.hpp:78
Teuchos_DefaultMpiComm.hpp
Panzer_ElementBlockIdToPhysicsIdMap.hpp
panzer::BlockedTpetraLinearObjFactory::getMap
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
Definition: Panzer_BlockedTpetraLinearObjFactory_impl.hpp:857
Panzer_IntrepidFieldPattern.hpp
Panzer_BlockedTpetraLinearObjFactory.hpp
panzer::UniqueGlobalIndexer::getElementBlockGIDCount
virtual int getElementBlockGIDCount(const std::string &blockId) const =0
How many GIDs are associate with a particular element block.
Panzer_WorksetContainer.hpp
panzer::L2Projection::buildRHSMatrix
Teuchos::RCP< Tpetra::CrsMatrix< double, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildRHSMatrix(const panzer::UniqueGlobalIndexer< LO, GO > &sourceDOFManager, const Teuchos::RCP< const Tpetra::Map< LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device >>> &ownedSourceMap, const std::string &sourceFieldName, const panzer::BasisDescriptor &sourceBasisDescriptor, const int vectorOrGradientDirectionIndex=-1)
Allocates, fills and returns a rectangular matrix for L2 projection of a scalar field, one dimension of gradient (for hgrad basis), or one dimension of a vector field onto the target scalar basis. If you wish to project all values of a vector field or all the gradients of a scalar field, then you will need three separate RHS matrices to form the RHS for each independently. The vectors must be independent Tpetra vectors to solve multiple right hand sides with the linear solver.
Definition: Panzer_L2Projection_impl.hpp:309
panzer::UniqueGlobalIndexer::getFieldNum
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
panzer::BlockedTpetraLinearObjFactory::getGhostedTpetraMatrix
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
Definition: Panzer_BlockedTpetraLinearObjFactory_impl.hpp:1064
Panzer_WorksetDescriptor.hpp
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
panzer::UniqueGlobalIndexer::getOwnedIndices
virtual void getOwnedIndices(std::vector< GlobalOrdinalT > &indices) const =0
Get the set of indices owned by this processor.
panzer::L2Projection::setup
void setup(const panzer::BasisDescriptor &targetBasis, const panzer::IntegrationDescriptor &integrationDescriptor, const Teuchos::RCP< const Teuchos::MpiComm< int >> &comm, const Teuchos::RCP< const panzer::ConnManager< LO, GO >> &connManager, const std::vector< std::string > &elementBlockNames, const Teuchos::RCP< panzer::WorksetContainer > worksetContainer=Teuchos::null)
Setup base objects for L2 Projections - requires target scalar basis and creates worksets if not supp...
Definition: Panzer_L2Projection_impl.hpp:26
panzer::L2Projection::buildMassMatrix
Teuchos::RCP< Tpetra::CrsMatrix< double, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildMassMatrix(bool use_lumping=false, const std::unordered_map< std::string, double > *elementBlockMultipliers=nullptr)
Allocates, fills and returns a mass matrix for L2 projection onto a target basis. ...
Definition: Panzer_L2Projection_impl.hpp:76
offsets
Kokkos::View< const int *, PHX::Device > offsets
Definition: Panzer_DOFCurl_impl.hpp:180
panzer::ALL_ELEMENTS
Workset size is set to the total number of local elements in the MPI process.
Definition: Panzer_WorksetDescriptor.hpp:58
Generated by   doxygen 1.8.5