Panzer Version of the Day
Loading...
Searching...
No Matches
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"
14#include "Panzer_TpetraLinearObjFactory.hpp"
22#include "Panzer_Workset.hpp"
23
24namespace 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->getLocalNumElements(),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));
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
PHX::View< const int * > offsets
Kokkos::View< typename Sacado::ScalarType< ScalarT >::type **, Kokkos::LayoutRight, PHX::Device > vals
const std::string & getType() const
Get type of basis.
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
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.
virtual void getOwnedAndGhostedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of owned and ghosted indices for this processor.
virtual int getElementBlockGIDCount(const std::size_t &blockIndex) const =0
How any GIDs are associate with each element in a particular element block.
virtual void getOwnedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of indices owned by this processor.
const Kokkos::View< const panzer::LocalOrdinal *, Kokkos::LayoutRight, PHX::Device > getElementLIDs(panzer::LocalOrdinal localElmtId) const
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
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,...
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.
Teuchos::RCP< panzer::GlobalIndexer > getTargetGlobalIndexer() const
Returns the target global indexer. Will be null if setup() has not been called.
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...
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...
@ ALL_ELEMENTS
Workset size is set to the total number of local elements in the MPI process.