Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_ScatterDirichletResidual_Epetra_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_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP
44#define PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP
45
46#include "Teuchos_RCP.hpp"
47#include "Teuchos_Assert.hpp"
48
49#include "Phalanx_DataLayout.hpp"
50
51#include "Epetra_Map.h"
52#include "Epetra_Vector.h"
53#include "Epetra_CrsMatrix.h"
54
56#include "Panzer_PureBasis.hpp"
61
62#include "Phalanx_DataLayout_MDALayout.hpp"
63
64#include "Teuchos_FancyOStream.hpp"
65
66// **********************************************************************
67// Specialization: Residual
68// **********************************************************************
69
70
71template<typename TRAITS,typename LO,typename GO>
73ScatterDirichletResidual_Epetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
74 const Teuchos::RCP<const panzer::GlobalIndexer> & /* cIndexer */,
75 const Teuchos::ParameterList& p)
76 : globalIndexer_(indexer)
77 , globalDataKey_("Residual Scatter Container")
78{
79 std::string scatterName = p.get<std::string>("Scatter Name");
80 scatterHolder_ =
81 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
82
83 // get names to be evaluated
84 const std::vector<std::string>& names =
85 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
86
87 // grab map from evaluated names to field names
88 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
89
90 // determine if we are scattering an initial condition
91 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
92
93 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
94 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
95 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
96 if (!scatterIC_) {
97 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
98 local_side_id_ = p.get<int>("Local Side ID");
99 }
100
101 // build the vector of fields that this is dependent on
102 scatterFields_.resize(names.size());
103 for (std::size_t eq = 0; eq < names.size(); ++eq) {
104 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
105
106 // tell the field manager that we depend on this field
107 this->addDependentField(scatterFields_[eq]);
108 }
109
110 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
111 if (checkApplyBC_) {
112 applyBC_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
114 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
115 this->addDependentField(applyBC_[eq]);
116 }
117 }
118
119 // this is what this evaluator provides
120 this->addEvaluatedField(*scatterHolder_);
121
122 if (p.isType<std::string>("Global Data Key"))
123 globalDataKey_ = p.get<std::string>("Global Data Key");
124
125 this->setName(scatterName+" Scatter Residual");
126}
127
128// **********************************************************************
129template<typename TRAITS,typename LO,typename GO>
131postRegistrationSetup(typename TRAITS::SetupData /* d */,
133{
134 fieldIds_.resize(scatterFields_.size());
135
136 // load required field numbers for fast use
137 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
138 // get field ID from DOF manager
139 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
140 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
141 }
142
143 // get the number of nodes (Should be renamed basis)
144 num_nodes = scatterFields_[0].extent(1);
145}
146
147// **********************************************************************
148template<typename TRAITS,typename LO,typename GO>
150preEvaluate(typename TRAITS::PreEvalData d)
151{
152 // extract linear object container
153 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
154
155 if(epetraContainer_==Teuchos::null) {
156 // extract linear object container
157 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
158 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
159
160 dirichletCounter_ = Teuchos::null;
161 }
162 else {
163 // extract dirichlet counter from container
164 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
165 = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject("Dirichlet Counter"),true);
166
167 dirichletCounter_ = epetraContainer->get_f();
168 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
169 }
170}
171
172// **********************************************************************
173template<typename TRAITS,typename LO,typename GO>
175evaluateFields(typename TRAITS::EvalData workset)
176{
177 // for convenience pull out some objects from workset
178 std::string blockId = this->wda(workset).block_id;
179 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
180
181 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
182 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
183 epetraContainer->get_f() :
184 epetraContainer->get_x();
185 // NOTE: A reordering of these loops will likely improve performance
186 // The "getGIDFieldOffsets may be expensive. However the
187 // "getElementGIDs" can be cheaper. However the lookup for LIDs
188 // may be more expensive!
189
190 // scatter operation for each cell in workset
191 auto LIDs = globalIndexer_->getLIDs();
192 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
193 Kokkos::deep_copy(LIDs_h, LIDs);
194
195
196 // loop over each field to be scattered
197 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
198 int fieldNum = fieldIds_[fieldIndex];
199 auto field = PHX::as_view(scatterFields_[fieldIndex]);
200 auto field_h = Kokkos::create_mirror_view(field);
201 Kokkos::deep_copy(field_h, field);
202
203 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
204 std::size_t cellLocalId = localCellIds[worksetCellIndex];
205
206 if (!scatterIC_) {
207 // this call "should" get the right ordering according to the Intrepid2 basis
208 const std::pair<std::vector<int>,std::vector<int> > & indicePair
209 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
210 const std::vector<int> & elmtOffset = indicePair.first;
211 const std::vector<int> & basisIdMap = indicePair.second;
212
213 // loop over basis functions
214 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
215 int offset = elmtOffset[basis];
216 int lid = LIDs_h(cellLocalId, offset);
217 if(lid<0) // not on this processor!
218 continue;
219
220 int basisId = basisIdMap[basis];
221
222 if (checkApplyBC_)
223 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
224 continue;
225
226 (*r)[lid] = field_h(worksetCellIndex,basisId);
227
228 // record that you set a dirichlet condition
229 if(dirichletCounter_!=Teuchos::null)
230 (*dirichletCounter_)[lid] = 1.0;
231 }
232 } else {
233 // this call "should" get the right ordering according to the Intrepid2 basis
234 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
235
236 // loop over basis functions
237 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
238 int offset = elmtOffset[basis];
239 int lid = LIDs_h(cellLocalId, offset);
240 if(lid<0) // not on this processor!
241 continue;
242
243 (*r)[lid] = field_h(worksetCellIndex,basis);
244
245 // record that you set a dirichlet condition
246 if(dirichletCounter_!=Teuchos::null)
247 (*dirichletCounter_)[lid] = 1.0;
248 }
249 }
250 }
251 }
252}
253
254// **********************************************************************
255// Specialization: Tangent
256// **********************************************************************
257
258
259template<typename TRAITS,typename LO,typename GO>
261ScatterDirichletResidual_Epetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
262 const Teuchos::RCP<const panzer::GlobalIndexer> & /* cIndexer */,
263 const Teuchos::ParameterList& p)
264 : globalIndexer_(indexer)
265 , globalDataKey_("Residual Scatter Container")
266{
267 std::string scatterName = p.get<std::string>("Scatter Name");
268 scatterHolder_ =
269 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
270
271 // get names to be evaluated
272 const std::vector<std::string>& names =
273 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
274
275 // grab map from evaluated names to field names
276 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
277
278 // determine if we are scattering an initial condition
279 scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
280
281 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
282 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
283 p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
284 if (!scatterIC_) {
285 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
286 local_side_id_ = p.get<int>("Local Side ID");
287 }
288
289 // build the vector of fields that this is dependent on
290 scatterFields_.resize(names.size());
291 for (std::size_t eq = 0; eq < names.size(); ++eq) {
292 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
293
294 // tell the field manager that we depend on this field
295 this->addDependentField(scatterFields_[eq]);
296 }
297
298 checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
299 if (checkApplyBC_) {
300 applyBC_.resize(names.size());
301 for (std::size_t eq = 0; eq < names.size(); ++eq) {
302 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
303 this->addDependentField(applyBC_[eq]);
304 }
305 }
306
307 // this is what this evaluator provides
308 this->addEvaluatedField(*scatterHolder_);
309
310 if (p.isType<std::string>("Global Data Key"))
311 globalDataKey_ = p.get<std::string>("Global Data Key");
312
313 this->setName(scatterName+" Scatter Tangent");
314}
315
316// **********************************************************************
317template<typename TRAITS,typename LO,typename GO>
319postRegistrationSetup(typename TRAITS::SetupData /* d */,
321{
322 fieldIds_.resize(scatterFields_.size());
323
324 // load required field numbers for fast use
325 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
326 // get field ID from DOF manager
327 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
328 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
329 }
330
331 // get the number of nodes (Should be renamed basis)
332 num_nodes = scatterFields_[0].extent(1);
333}
334
335// **********************************************************************
336template<typename TRAITS,typename LO,typename GO>
338preEvaluate(typename TRAITS::PreEvalData d)
339{
340 // extract linear object container
341 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
342
343 if(epetraContainer_==Teuchos::null) {
344 // extract linear object container
345 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
346 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
347
348 dirichletCounter_ = Teuchos::null;
349 }
350 else {
351 // extract dirichlet counter from container
352 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
353 = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject("Dirichlet Counter"),true);
354
355 dirichletCounter_ = epetraContainer->get_f();
356 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
357 }
358
359 using Teuchos::RCP;
360 using Teuchos::rcp_dynamic_cast;
361
362 // this is the list of parameters and their names that this scatter has to account for
363 std::vector<std::string> activeParameters =
364 rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
365
366 // ETP 02/03/16: This code needs to be updated to properly handle scatterIC_
367 TEUCHOS_ASSERT(!scatterIC_);
368 dfdp_vectors_.clear();
369 for(std::size_t i=0;i<activeParameters.size();i++) {
370 RCP<Epetra_Vector> vec = rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
371 dfdp_vectors_.push_back(vec);
372 }
373}
374
375// **********************************************************************
376template<typename TRAITS,typename LO,typename GO>
378evaluateFields(typename TRAITS::EvalData workset)
379{
380 // for convenience pull out some objects from workset
381 std::string blockId = this->wda(workset).block_id;
382 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
383
384 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
385 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
386 epetraContainer->get_f() :
387 epetraContainer->get_x();
388 // NOTE: A reordering of these loops will likely improve performance
389 // The "getGIDFieldOffsets may be expensive. However the
390 // "getElementGIDs" can be cheaper. However the lookup for LIDs
391 // may be more expensive!
392
393 // loop over each field to be scattered
394 auto LIDs = globalIndexer_->getLIDs();
395 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
396 Kokkos::deep_copy(LIDs_h, LIDs);
397 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
398 int fieldNum = fieldIds_[fieldIndex];
399 auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
400 Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
401
402 // scatter operation for each cell in workset
403 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
404 std::size_t cellLocalId = localCellIds[worksetCellIndex];
405
406 if (!scatterIC_) {
407 // this call "should" get the right ordering according to the Intrepid2 basis
408 const std::pair<std::vector<int>,std::vector<int> > & indicePair
409 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
410 const std::vector<int> & elmtOffset = indicePair.first;
411 const std::vector<int> & basisIdMap = indicePair.second;
412
413 // loop over basis functions
414 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
415 int offset = elmtOffset[basis];
416 int lid = LIDs_h(cellLocalId, offset);
417 if(lid<0) // not on this processor!
418 continue;
419
420 int basisId = basisIdMap[basis];
421
422 if (checkApplyBC_)
423 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
424 continue;
425
426 ScalarT value = scatterField_h(worksetCellIndex,basisId);
427 // (*r)[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
428
429 // then scatter the sensitivity vectors
430 if(value.size()==0)
431 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
432 (*dfdp_vectors_[d])[lid] = 0.0;
433 else
434 for(int d=0;d<value.size();d++) {
435 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
436 }
437
438 // record that you set a dirichlet condition
439 if(dirichletCounter_!=Teuchos::null)
440 (*dirichletCounter_)[lid] = 1.0;
441 }
442 } else {
443 // this call "should" get the right ordering according to the Intrepid2 basis
444 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
445
446 // loop over basis functions
447 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
448 int offset = elmtOffset[basis];
449 int lid = LIDs_h(cellLocalId, offset);
450 if(lid<0) // not on this processor!
451 continue;
452
453 ScalarT value = scatterField_h(worksetCellIndex,basis);
454 // (*r)[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
455
456 // then scatter the sensitivity vectors
457 if(value.size()==0)
458 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
459 (*dfdp_vectors_[d])[lid] = 0.0;
460 else
461 for(int d=0;d<value.size();d++) {
462 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
463 }
464
465 // record that you set a dirichlet condition
466 if(dirichletCounter_!=Teuchos::null)
467 (*dirichletCounter_)[lid] = 1.0;
468 }
469 }
470 }
471 }
472}
473
474// **********************************************************************
475// Specialization: Jacobian
476// **********************************************************************
477
478template<typename TRAITS,typename LO,typename GO>
480ScatterDirichletResidual_Epetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
481 const Teuchos::RCP<const panzer::GlobalIndexer> & cIndexer,
482 const Teuchos::ParameterList& p)
483 : globalIndexer_(indexer)
484 , colGlobalIndexer_(cIndexer)
485 , globalDataKey_("Residual Scatter Container")
486 , preserveDiagonal_(false)
487{
488 std::string scatterName = p.get<std::string>("Scatter Name");
489 scatterHolder_ =
490 Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
491
492 // get names to be evaluated
493 const std::vector<std::string>& names =
494 *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
495
496 // grab map from evaluated names to field names
497 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
498
499 Teuchos::RCP<PHX::DataLayout> dl =
500 p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
501
502 side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
503 local_side_id_ = p.get<int>("Local Side ID");
504
505 // build the vector of fields that this is dependent on
506 scatterFields_.resize(names.size());
507 for (std::size_t eq = 0; eq < names.size(); ++eq) {
508 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
509
510 // tell the field manager that we depend on this field
511 this->addDependentField(scatterFields_[eq]);
512 }
513
514 checkApplyBC_ = p.get<bool>("Check Apply BC");
515 if (checkApplyBC_) {
516 applyBC_.resize(names.size());
517 for (std::size_t eq = 0; eq < names.size(); ++eq) {
518 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
519 this->addDependentField(applyBC_[eq]);
520 }
521 }
522
523 // this is what this evaluator provides
524 this->addEvaluatedField(*scatterHolder_);
525
526 if (p.isType<std::string>("Global Data Key"))
527 globalDataKey_ = p.get<std::string>("Global Data Key");
528
529 if (p.isType<bool>("Preserve Diagonal"))
530 preserveDiagonal_ = p.get<bool>("Preserve Diagonal");
531
532 this->setName(scatterName+" Scatter Residual (Jacobian)");
533}
534
535// **********************************************************************
536template<typename TRAITS,typename LO,typename GO>
538postRegistrationSetup(typename TRAITS::SetupData /* d */,
540{
541 fieldIds_.resize(scatterFields_.size());
542
543 // load required field numbers for fast use
544 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
545 // get field ID from DOF manager
546 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
547 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
548 }
549
550 // get the number of nodes (Should be renamed basis)
551 num_nodes = scatterFields_[0].extent(1);
552 num_eq = scatterFields_.size();
553}
554
555// **********************************************************************
556template<typename TRAITS,typename LO,typename GO>
558preEvaluate(typename TRAITS::PreEvalData d)
559{
560 // extract linear object container
561 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
562
563 if(epetraContainer_==Teuchos::null) {
564 // extract linear object container
565 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
566 epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc,true);
567
568 dirichletCounter_ = Teuchos::null;
569 }
570 else {
571 // extract dirichlet counter from container
572 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc->getDataObject("Dirichlet Counter");
573 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(dataContainer,true);
574
575 dirichletCounter_ = epetraContainer->get_f();
576 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
577 }
578}
579
580// **********************************************************************
581template<typename TRAITS,typename LO,typename GO>
583evaluateFields(typename TRAITS::EvalData workset)
584{
585 Kokkos::View<const int*, Kokkos::LayoutRight, PHX::Device> cLIDs, rLIDs;
586 int gidCount(0);
587 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
588 const Teuchos::RCP<const panzer::GlobalIndexer>&
589 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
590
591 // for convenience pull out some objects from workset
592 std::string blockId = this->wda(workset).block_id;
593 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
594
595 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
596 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
597 Teuchos::RCP<Epetra_Vector> r = epetraContainer->get_f();
598 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
599 // NOTE: A reordering of these loops will likely improve performance
600 // The "getGIDFieldOffsets may be expensive. However the
601 // "getElementGIDs" can be cheaper. However the lookup for LIDs
602 // may be more expensive!
603
604 // scatter operation for each cell in workset
605
606 auto LIDs = globalIndexer_->getLIDs();
607 auto colLIDs = colGlobalIndexer->getLIDs();
608 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
609 auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
610 Kokkos::deep_copy(LIDs_h, LIDs);
611 Kokkos::deep_copy(colLIDs_h, colLIDs);
612
613 std::vector<typename decltype(scatterFields_[0].get_static_view())::HostMirror> scatterFields_h;
614 for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
615 scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
616 Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
617 }
618
619 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
620 std::size_t cellLocalId = localCellIds[worksetCellIndex];
621
622 gidCount = colGlobalIndexer->getElementBlockGIDCount(blockId);
623
624 // loop over each field to be scattered
625 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
626 int fieldNum = fieldIds_[fieldIndex];
627
628 // this call "should" get the right ordering according to the Intrepid2 basis
629 const std::pair<std::vector<int>,std::vector<int> > & indicePair
630 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
631 const std::vector<int> & elmtOffset = indicePair.first;
632 const std::vector<int> & basisIdMap = indicePair.second;
633
634 // loop over basis functions
635 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
636 int offset = elmtOffset[basis];
637 int row = LIDs_h(cellLocalId, offset);
638 if(row<0) // not on this processor
639 continue;
640
641 int basisId = basisIdMap[basis];
642
643 if (checkApplyBC_)
644 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
645 continue;
646
647 // zero out matrix row
648 {
649 int numEntries = 0;
650 int * rowIndices = 0;
651 double * rowValues = 0;
652
653 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
654
655 for(int i=0;i<numEntries;i++) {
656 if(preserveDiagonal_) {
657 if(row!=rowIndices[i])
658 rowValues[i] = 0.0;
659 }
660 else
661 rowValues[i] = 0.0;
662 }
663 }
664
665 // int gid = GIDs[offset];
666 const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,basisId);
667
668 if(r!=Teuchos::null)
669 (*r)[row] = scatterField.val();
670 if(dirichletCounter_!=Teuchos::null) {
671 // std::cout << "Writing " << row << " " << dirichletCounter_->MyLength() << std::endl;
672 (*dirichletCounter_)[row] = 1.0; // mark row as dirichlet
673 }
674
675 // loop over the sensitivity indices: all DOFs on a cell
676 std::vector<double> jacRow(scatterField.size(),0.0);
677
678 if(!preserveDiagonal_) {
679 int err = Jac->ReplaceMyValues(row, gidCount, scatterField.dx(),
680 &colLIDs_h(cellLocalId,0));
681 TEUCHOS_ASSERT(err==0);
682 }
683 }
684 }
685 }
686}
687
688// **********************************************************************
689
690#endif
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
Pushes residual values into the residual vector for a Newton-based solve.