Panzer  Version of the Day
Panzer_L2Projection.cpp
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_L2Projection.hpp"
14 #include "Panzer_TpetraLinearObjFactory.hpp"
22 #include "Panzer_Workset.hpp"
23 
24 namespace panzer {
25 
27  const panzer::IntegrationDescriptor& integrationDescriptor,
28  const Teuchos::RCP<const Teuchos::MpiComm<int>>& comm,
29  const Teuchos::RCP<const panzer::ConnManager>& 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(Teuchos::rcp_const_cast<panzer::ConnManager>(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  Teuchos::RCP<panzer::GlobalIndexer>
71  {return targetGlobalIndexer_;}
72 
73  Teuchos::RCP<Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>
75  const std::unordered_map<std::string,double>* elementBlockMultipliers)
76  {
77  PANZER_FUNC_TIME_MONITOR_DIFF("L2Projection::Build Mass Matrix",TopBuildMassMatrix);
78  TEUCHOS_ASSERT(setupCalled_);
79 
80  if (elementBlockMultipliers != nullptr) {
81  TEUCHOS_ASSERT(elementBlockMultipliers->size() == elementBlockNames_.size());
82  }
83 
84  // Allocate the owned matrix
85  std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> indexers;
86  indexers.push_back(targetGlobalIndexer_);
87 
89 
90  auto ownedMatrix = factory.getTpetraMatrix(0,0);
91  auto ghostedMatrix = factory.getGhostedTpetraMatrix(0,0);
92  ownedMatrix->resumeFill();
93  ownedMatrix->setAllToScalar(0.0);
94  ghostedMatrix->resumeFill();
95  ghostedMatrix->setAllToScalar(0.0);
96 
97  const int fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
98 
99  const bool is_scalar = targetBasisDescriptor_.getType()=="HGrad" || targetBasisDescriptor_.getType()=="Const" || targetBasisDescriptor_.getType()=="HVol";
100 
101  // Loop over element blocks and fill mass matrix
102  if(is_scalar){
103  auto M = ghostedMatrix->getLocalMatrixDevice();
104  for (const auto& block : elementBlockNames_) {
105 
106  double ebMultiplier = 1.0;
107  if (elementBlockMultipliers != nullptr)
108  ebMultiplier = elementBlockMultipliers->find(block)->second;
109 
110  // Based on descriptor, currently assumes there should only be one workset
112  const auto worksets = worksetContainer_->getWorksets(wd);
113 
114  for (const auto& workset : *worksets) {
115 
116  const auto basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
117 
118  const auto unweightedBasis = basisValues.getBasisValues(false).get_static_view();
119  const auto weightedBasis = basisValues.getBasisValues(true).get_static_view();
120 
121  const std::vector<panzer::LocalOrdinal>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
122  PHX::View<panzer::LocalOrdinal*> kOffsets("MassMatrix: Offsets",offsets.size());
123  auto kOffsets_h = Kokkos::create_mirror_view(kOffsets);
124 
125  for(const auto& i : offsets)
126  kOffsets_h(i) = offsets[i];
127 
128  Kokkos::deep_copy(kOffsets, kOffsets_h);
129 
130  // Local Ids
131  PHX::View<panzer::LocalOrdinal**> localIds("MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
132  targetGlobalIndexer_->getElementBlockGIDCount(block));
133 
134  // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
135  const auto cellLocalIdsNoGhost = Kokkos::subview(workset.getLocalCellIDs(),std::make_pair(0,workset.numOwnedCells()));
136 
137  targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
138 
139  const int numBasisPoints = static_cast<int>(weightedBasis.extent(1));
140  if ( use_lumping ) {
141  Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (const int& cell) {
142  double total_mass = 0.0, trace = 0.0;
143 
144  panzer::LocalOrdinal cLIDs[256];
145  const int numIds = static_cast<int>(localIds.extent(1));
146  for(int i=0;i<numIds;++i)
147  cLIDs[i] = localIds(cell,i);
148 
149  double vals[256]={0.0};
150  const int numQP = static_cast<int>(unweightedBasis.extent(2));
151 
152  for (int row=0; row < numBasisPoints; ++row) {
153  for (int col=0; col < numIds; ++col) {
154  for (int qp=0; qp < numQP; ++qp) {
155  auto tmp = unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
156  total_mass += tmp;
157  if (col == row )
158  trace += tmp;
159  }
160  }
161  }
162 
163  for (int row=0; row < numBasisPoints; ++row) {
164  for (int col=0; col < numBasisPoints; ++col)
165  vals[col] = 0;
166 
167  int offset = kOffsets(row);
168  panzer::LocalOrdinal lid = localIds(cell,offset);
169  int col = row;
170  vals[col] = 0.0;
171  for (int qp=0; qp < numQP; ++qp)
172  vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier * total_mass / trace;
173 
174  M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
175  }
176  });
177 
178  } else {
179  Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (const int& cell) {
180 
181  panzer::LocalOrdinal cLIDs[256];
182  const int numIds = static_cast<int>(localIds.extent(1));
183  for(int i=0;i<numIds;++i)
184  cLIDs[i] = localIds(cell,i);
185 
186  double vals[256];
187  const int numQP = static_cast<int>(unweightedBasis.extent(2));
188 
189  for (int row=0; row < numBasisPoints; ++row) {
190  int offset = kOffsets(row);
191  panzer::LocalOrdinal lid = localIds(cell,offset);
192 
193  for (int col=0; col < numIds; ++col) {
194  vals[col] = 0.0;
195  for (int qp=0; qp < numQP; ++qp)
196  vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
197  }
198  M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
199 
200  }
201 
202  });
203  }
204  }
205  }
206  } else {
207  auto M = ghostedMatrix->getLocalMatrixDevice();
208  for (const auto& block : elementBlockNames_) {
209 
210  double ebMultiplier = 1.0;
211  if (elementBlockMultipliers != nullptr)
212  ebMultiplier = elementBlockMultipliers->find(block)->second;
213 
214  // Based on descriptor, currently assumes there should only be one workset
216  const auto& worksets = worksetContainer_->getWorksets(wd);
217 
218  for (const auto& workset : *worksets) {
219 
220  const auto basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
221 
222  const auto unweightedBasis = basisValues.getVectorBasisValues(false).get_static_view();
223  const auto weightedBasis = basisValues.getVectorBasisValues(true).get_static_view();
224 
225  const std::vector<panzer::LocalOrdinal>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
226  PHX::View<panzer::LocalOrdinal*> kOffsets("MassMatrix: Offsets",offsets.size());
227  auto kOffsets_h = Kokkos::create_mirror_view(kOffsets);
228 
229  for(const auto& i : offsets)
230  kOffsets_h(i) = offsets[i];
231 
232  Kokkos::deep_copy(kOffsets, kOffsets_h);
233 
234  // Local Ids
235  PHX::View<panzer::LocalOrdinal**> 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 PHX::View<const int*> 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  panzer::LocalOrdinal 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  panzer::LocalOrdinal 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 
274  {
275  PANZER_FUNC_TIME_MONITOR_DIFF("Exporting of mass matrix",ExportMM);
276  auto map = factory.getMap(0);
277  ghostedMatrix->fillComplete(map,map);
278  const auto exporter = factory.getGhostedExport(0);
279  ownedMatrix->doExport(*ghostedMatrix, *exporter, Tpetra::ADD);
280  ownedMatrix->fillComplete();
281  }
282  return ownedMatrix;
283  }
284 
285  Teuchos::RCP<Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>
287  {
288  PANZER_FUNC_TIME_MONITOR_DIFF("L2Projection<panzer::LocalOrdinal,panzer::GlobalOrdinal>::buildInverseLumpedMassMatrix",buildInvLMM);
289  using Teuchos::rcp;
290  const auto massMatrix = this->buildMassMatrix(true);
291  const auto lumpedMassMatrix = rcp(new Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>(massMatrix->getDomainMap(),1,true));
292  const auto tmp = rcp(new Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>(massMatrix->getRangeMap(),1,false));
293  tmp->putScalar(1.0);
294  {
295  PANZER_FUNC_TIME_MONITOR_DIFF("Apply",Apply);
296  massMatrix->apply(*tmp,*lumpedMassMatrix);
297  }
298  {
299  PANZER_FUNC_TIME_MONITOR_DIFF("Reciprocal",Reciprocal);
300  lumpedMassMatrix->reciprocal(*lumpedMassMatrix);
301  }
302  return lumpedMassMatrix;
303  }
304 
305  Teuchos::RCP<Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>
307  const Teuchos::RCP<const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>& inputOwnedSourceMap,
308  const std::string& sourceFieldName,
309  const panzer::BasisDescriptor& sourceBasisDescriptor,
310  const int directionIndex)
311  {
312  // *******************
313  // Create Retangular matrix (both ghosted and owned).
314  // *******************
315  using Teuchos::RCP;
316  using Teuchos::rcp;
317  using MapType = Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
318  using GraphType = Tpetra::CrsGraph<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
319  using ExportType = Tpetra::Export<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
320  using MatrixType = Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
321 
322  // *******************
323  // Build ghosted graph
324  // *******************
325 
326  RCP<MapType> ghostedTargetMap;
327  {
328  std::vector<panzer::GlobalOrdinal> indices;
329  targetGlobalIndexer_->getOwnedAndGhostedIndices(indices);
330  ghostedTargetMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
331  }
332 
333  RCP<MapType> ghostedSourceMap;
334  {
335  std::vector<panzer::GlobalOrdinal> indices;
336  sourceGlobalIndexer.getOwnedAndGhostedIndices(indices);
337  ghostedSourceMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
338  }
339 
340 
341  // Now insert the non-zero pattern per row
342  // count number of entries per row; required by CrsGraph constructor
343  std::vector<size_t> nEntriesPerRow(ghostedTargetMap->getNodeNumElements(),0);
344  std::vector<std::string> elementBlockIds;
345  targetGlobalIndexer_->getElementBlockIds(elementBlockIds);
346  std::vector<std::string>::const_iterator blockItr;
347  for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
348  std::string blockId = *blockItr;
349  const std::vector<panzer::LocalOrdinal> & elements = targetGlobalIndexer_->getElementBlock(blockId);
350 
351  std::vector<panzer::GlobalOrdinal> row_gids;
352  std::vector<panzer::GlobalOrdinal> col_gids;
353 
354  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
355  targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
356  sourceGlobalIndexer.getElementGIDs(elements[elmt],col_gids);
357  for(std::size_t row=0;row<row_gids.size();row++) {
358  panzer::LocalOrdinal lid =
359  ghostedTargetMap->getLocalElement(row_gids[row]);
360  nEntriesPerRow[lid] += col_gids.size();
361  }
362  }
363  }
364 
365  Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
366  RCP<GraphType> ghostedGraph = rcp(new GraphType(ghostedTargetMap,ghostedSourceMap,nEntriesPerRowView,Tpetra::StaticProfile));
367 
368  for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
369  std::string blockId = *blockItr;
370  const std::vector<panzer::LocalOrdinal> & elements = targetGlobalIndexer_->getElementBlock(blockId);
371 
372  std::vector<panzer::GlobalOrdinal> row_gids;
373  std::vector<panzer::GlobalOrdinal> col_gids;
374 
375  for(std::size_t elmt=0;elmt<elements.size();elmt++) {
376  targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
377  sourceGlobalIndexer.getElementGIDs(elements[elmt],col_gids);
378  for(std::size_t row=0;row<row_gids.size();row++)
379  ghostedGraph->insertGlobalIndices(row_gids[row],col_gids);
380  }
381  }
382 
383  RCP<MapType> ownedTargetMap;
384  {
385  std::vector<panzer::GlobalOrdinal> indices;
386  targetGlobalIndexer_->getOwnedIndices(indices);
387  ownedTargetMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
388  }
389 
390  RCP<const MapType> ownedSourceMap = inputOwnedSourceMap;
391  if (is_null(ownedSourceMap)) {
392  std::vector<panzer::GlobalOrdinal> indices;
393  sourceGlobalIndexer.getOwnedIndices(indices);
394  ownedSourceMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
395  }
396 
397  // Fill complete with owned range and domain map
398  ghostedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
399 
400  // *****************
401  // Build owned graph
402  // *****************
403 
404  RCP<GraphType> ownedGraph = rcp(new GraphType(ownedTargetMap,0));
405  RCP<const ExportType> exporter = rcp(new ExportType(ghostedTargetMap,ownedTargetMap));
406  ownedGraph->doExport(*ghostedGraph, *exporter, Tpetra::INSERT);
407  ownedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
408 
409  RCP<MatrixType> ghostedMatrix = rcp(new MatrixType(ghostedGraph));
410  RCP<MatrixType> ownedMatrix = rcp(new MatrixType(ownedGraph));
411  // ghostedMatrix.fillComplete();
412  // ghostedMatrix.resumeFill();
413 
414  ghostedMatrix->setAllToScalar(0.0);
415  ownedMatrix->setAllToScalar(0.0);
416 
417  // *******************
418  // Fill ghosted matrix
419  // *******************
420  for (const auto& block : elementBlockNames_) {
421 
423  const auto& worksets = worksetContainer_->getWorksets(wd);
424  for (const auto& workset : *worksets) {
425 
426  // Get target basis values: current implementation assumes target basis is HGrad
427  const auto& targetBasisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
428  const auto& targetWeightedBasis = targetBasisValues.getBasisValues(true).get_static_view();
429 
430  // Sources can be any basis
431  const auto& sourceBasisValues = workset.getBasisValues(sourceBasisDescriptor,integrationDescriptor_);
432  PHX::View<const double***> sourceUnweightedScalarBasis;
433  PHX::View<const double****> sourceUnweightedVectorBasis;
434  bool useRankThreeBasis = false; // default to gradient or vector basis
435  if ( (sourceBasisDescriptor.getType() == "HGrad") || (sourceBasisDescriptor.getType() == "Const") || (sourceBasisDescriptor.getType() == "HVol") ) {
436  if (directionIndex == -1) { // Project dof value
437  sourceUnweightedScalarBasis = sourceBasisValues.getBasisValues(false).get_static_view();
438  useRankThreeBasis = true;
439  }
440  else { // Project dof gradient of scalar basis
441  sourceUnweightedVectorBasis = sourceBasisValues.getGradBasisValues(false).get_static_view();
442  }
443  }
444  else { // Project vector value
445  sourceUnweightedVectorBasis = sourceBasisValues.getVectorBasisValues(false).get_static_view();
446  }
447 
448  // Get the element local ids
449  PHX::View<panzer::LocalOrdinal**> targetLocalIds("buildRHSMatrix: targetLocalIds", workset.numOwnedCells(),
450  targetGlobalIndexer_->getElementBlockGIDCount(block));
451  PHX::View<panzer::LocalOrdinal**> sourceLocalIds("buildRHSMatrix: sourceLocalIds", workset.numOwnedCells(),
452  sourceGlobalIndexer.getElementBlockGIDCount(block));
453  {
454  // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
455  const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
456  targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,targetLocalIds);
457  sourceGlobalIndexer.getElementLIDs(cellLocalIdsNoGhost,sourceLocalIds);
458  }
459 
460  // Get the offsets
461  PHX::View<panzer::LocalOrdinal*> targetFieldOffsets;
462  {
463  const auto fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
464  const std::vector<panzer::LocalOrdinal>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
465  targetFieldOffsets = PHX::View<panzer::LocalOrdinal*>("L2Projection:buildRHS:targetFieldOffsets",offsets.size());
466  const auto hostOffsets = Kokkos::create_mirror_view(targetFieldOffsets);
467  for(size_t i=0; i < offsets.size(); ++i)
468  hostOffsets(i) = offsets[i];
469  Kokkos::deep_copy(targetFieldOffsets,hostOffsets);
470  }
471 
472  PHX::View<panzer::LocalOrdinal*> sourceFieldOffsets;
473  {
474  const auto fieldIndex = sourceGlobalIndexer.getFieldNum(sourceFieldName);
475  const std::vector<panzer::LocalOrdinal>& offsets = sourceGlobalIndexer.getGIDFieldOffsets(block,fieldIndex);
476  TEUCHOS_TEST_FOR_EXCEPTION(offsets.size() == 0, std::runtime_error,
477  "ERROR: panzer::L2Projection::buildRHSMatrix() - The source field, \""
478  << sourceFieldName << "\", does not exist in element block \""
479  << block << "\"!");
480  sourceFieldOffsets = PHX::View<panzer::LocalOrdinal*>("L2Projection:buildRHS:sourceFieldOffsets",offsets.size());
481  const auto hostOffsets = Kokkos::create_mirror_view(sourceFieldOffsets);
482  for(size_t i=0; i <offsets.size(); ++i)
483  hostOffsets(i) = offsets[i];
484  Kokkos::deep_copy(sourceFieldOffsets,hostOffsets);
485  }
486 
487  const auto localMatrix = ghostedMatrix->getLocalMatrixDevice();
488  const int numRows = static_cast<int>(targetWeightedBasis.extent(1));
489  int tmpNumCols = -1;
490  int tmpNumQP = -1;
491  if (useRankThreeBasis) {
492  tmpNumCols = static_cast<int>(sourceUnweightedScalarBasis.extent(1));
493  tmpNumQP = static_cast<int>(sourceUnweightedScalarBasis.extent(2));
494  }
495  else {
496  tmpNumCols = static_cast<int>(sourceUnweightedVectorBasis.extent(1));
497  tmpNumQP = static_cast<int>(sourceUnweightedVectorBasis.extent(2));
498  }
499  const int numCols = tmpNumCols;
500  const int numQP = tmpNumQP;
501 
502  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.numOwnedCells()),KOKKOS_LAMBDA (const int& cell) {
503  panzer::LocalOrdinal cLIDs[256];
504  double vals[256];
505  for (int row = 0; row < numRows; ++row) {
506  const int rowOffset = targetFieldOffsets(row);
507  const int rowLID = targetLocalIds(cell,rowOffset);
508  for (int col = 0; col < numCols; ++col)
509  vals[col] = 0.0;
510 
511  for (int col = 0; col < numCols; ++col) {
512  for (int qp = 0; qp < numQP; ++qp) {
513  const int colOffset = sourceFieldOffsets(col);
514  const int colLID = sourceLocalIds(cell,colOffset);
515  cLIDs[col] = colLID;
516  if (useRankThreeBasis)
517  vals[col] += sourceUnweightedScalarBasis(cell,col,qp) * targetWeightedBasis(cell,row,qp);
518  else
519  vals[col] += sourceUnweightedVectorBasis(cell,col,qp,directionIndex) * targetWeightedBasis(cell,row,qp);
520  }
521  }
522  localMatrix.sumIntoValues(rowLID,cLIDs,numCols,vals,true,true);
523  }
524  });
525  }
526 
527  }
528  ghostedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
529  ownedMatrix->resumeFill();
530  ownedMatrix->doExport(*ghostedMatrix,*exporter,Tpetra::ADD);
531  ownedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
532 
533  return ownedMatrix;
534  }
535 
536 }
537 
538 #endif
void setup(const panzer::BasisDescriptor &targetBasis, const panzer::IntegrationDescriptor &integrationDescriptor, const Teuchos::RCP< const Teuchos::MpiComm< int >> &comm, const Teuchos::RCP< const panzer::ConnManager > &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...
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
virtual int getElementBlockGIDCount(const std::size_t &blockIndex) const =0
How any GIDs are associate with a particular element block.
Teuchos::RCP< panzer::GlobalIndexer > getTargetGlobalIndexer() const
Returns the target global indexer. Will be null if setup() has not been called.
PHX::View< const int * > offsets
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
virtual void getOwnedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of indices owned by this processor.
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
virtual void getOwnedAndGhostedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of owned and ghosted indices for this processor.
Teuchos::RCP< Tpetra::CrsMatrix< double, panzer::LocalOrdinal, panzer::GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildRHSMatrix(const panzer::GlobalIndexer &sourceDOFManager, const Teuchos::RCP< const Tpetra::Map< panzer::LocalOrdinal, panzer::GlobalOrdinal, 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.
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.
virtual void getElementGIDs(panzer::LocalOrdinal localElmtId, std::vector< panzer::GlobalOrdinal > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
Teuchos::RCP< Tpetra::MultiVector< double, panzer::LocalOrdinal, panzer::GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildInverseLumpedMassMatrix()
Allocates, fills and returns a Tpetra::MultiVector containing the inverse lumped mass matrix values a...
const std::string & getType() const
Get type of basis.
const Kokkos::View< const panzer::LocalOrdinal *, Kokkos::LayoutRight, PHX::Device > getElementLIDs(panzer::LocalOrdinal localElmtId) const
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
Workset size is set to the total number of local elements in the MPI process.
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
Teuchos::RCP< Tpetra::CrsMatrix< double, panzer::LocalOrdinal, panzer::GlobalOrdinal, 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. ...