1#ifndef BCBuilder_def_hpp
2#define BCBuilder_def_hpp
4#include "feddlib/core/Utils/FEDDUtils.hpp"
5#include <Teuchos_RCPDecl.hpp>
6#include <Teuchos_ScalarTraitsDecl.hpp>
9inline void dummyFuncBC(
double* x,
double* res,
double t,
const double* parameters)
24template<
class SC,
class LO,
class GO,
class NO>
25BCBuilder<SC,LO,GO,NO>::BCBuilder():
37,SetSystemRowTimer_(Teuchos::TimeMonitor::getNewCounter(
"BCBuilder: SetSystemRow")),
38BlockRowHasDirichletTimer_(Teuchos::TimeMonitor::getNewCounter(
"BCBuilder: SetSystemRow: BlockHasDiri")),
39FindFlagTimer_(Teuchos::TimeMonitor::getNewCounter(
"BCBuilder: SetSystemRow: BlockHasDiri: FindFlag"))
45template<
class SC,
class LO,
class GO,
class NO>
48 vecBC_func_.push_back(funcBC);
49 vecFlag_.push_back(flag);
50 vecBlockID_.push_back(block);
51 vecDomain_.push_back(domain);
52 vecBCType_.push_back(type);
53 vecDofs_.push_back(dofs);
54 vec_dbl_Type dummy(1,0.);
55 vecBC_Parameters_.push_back(dummy);
56 MultiVectorConstPtr_Type dummyMV;
57 vecExternalSol_.push_back(dummyMV);
58 vecFlowRateBool_.push_back(
false);
59 vecBC_func_flowRate_.push_back(dummyFuncBC);
64template<
class SC,
class LO,
class GO,
class NO>
65void BCBuilder<SC,LO,GO,NO>::addBC(BC_func_Type funcBC,
int flag,
int block,
const DomainPtr_Type &domain, std::string type,
int dofs, vec_dbl_Type ¶meter_vec){
67 vecBC_func_.push_back(funcBC);
68 vecFlag_.push_back(flag);
69 vecBlockID_.push_back(block);
70 vecDomain_.push_back(domain);
71 vecBCType_.push_back(type);
72 vecDofs_.push_back(dofs);
73 vecBC_Parameters_.push_back(parameter_vec);
74 MultiVectorConstPtr_Type dummyMV;
75 vecExternalSol_.push_back(dummyMV);
76 vecFlowRateBool_.push_back(
false);
77 vecBC_func_flowRate_.push_back(dummyFuncBC);
81template<
class SC,
class LO,
class GO,
class NO>
82void BCBuilder<SC,LO,GO,NO>::addBC(BC_func_Type funcBC,
int flag,
int block,
const DomainPtr_Type &domain, std::string type,
int dofs, vec_dbl_Type ¶meter_vec, MultiVectorConstPtr_Type& externalSol){
84 vecBC_func_.push_back(funcBC);
85 vecFlag_.push_back(flag);
86 vecBlockID_.push_back(block);
87 vecDomain_.push_back(domain);
88 vecBCType_.push_back(type);
89 vecDofs_.push_back(dofs);
90 vecBC_Parameters_.push_back(parameter_vec);
91 vecExternalSol_.push_back(externalSol);
92 vecFlowRateBool_.push_back(
false);
93 vecBC_func_flowRate_.push_back(dummyFuncBC);
96template<
class SC,
class LO,
class GO,
class NO>
98 std::string type,
int dofs, vec_dbl_Type ¶meter_vec, MultiVectorConstPtr_Type& externalSol,
99 bool determineFlowRate,BC_func_Type funcBC_flowRate){
101 vecBC_func_.push_back(funcBC);
102 vecFlag_.push_back(flag);
103 vecBlockID_.push_back(block);
104 vecDomain_.push_back(domain);
105 vecBCType_.push_back(type);
106 vecDofs_.push_back(dofs);
107 vecBC_Parameters_.push_back(parameter_vec);
108 vecExternalSol_.push_back(externalSol);
109 vecFlowRateBool_.push_back(
true);
110 vecBC_func_flowRate_.push_back(funcBC_flowRate);
113template<
class SC,
class LO,
class GO,
class NO>
122template<
class SC,
class LO,
class GO,
class NO>
127 TEUCHOS_TEST_FOR_EXCEPTION( blockMV->getNumVectors()>1, std::runtime_error,
"BCBuilder::setRHS() only for getNumVectors == 1.");
130 for (
int block = 0; block < blockMV->size(); block++) {
132 if (vecFlowRateBool_[loc]) {
133 this->determineVelocityForFlowrate(loc,t);
135 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
136 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
137 int dim = vecDomain_.at(loc)->getDimension();
138 int dofs = vecDofs_.at(loc);
140 result.resize(dofs,0.);
141 point.resize(dim,0.);
143 resultPtr_ = Teuchos::rcp(
new vec_dbl_Type(dofs,0.));
144 pointPtr_ = Teuchos::rcp(
new vec_dbl_Type(dim,0.));
146 Teuchos::ArrayRCP<SC> valuesRHS = blockMV->getBlock(block)->getDataNonConst(0);
147 for (
int i = 0 ; i < bcFlags->size(); i++) {
148 int flag = bcFlags->at(i);
150 if (!vecBCType_.at(loc).compare(
"Dirichlet") ||
151 !vecBCType_.at(loc).compare(
"Dirichlet_X") ||
152 !vecBCType_.at(loc).compare(
"Dirichlet_Y") ||
153 !vecBCType_.at(loc).compare(
"Dirichlet_Z")||
154 !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ||
155 !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ||
156 !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z")) {
157 if (!vecExternalSol_[loc].is_null()) {
158 std::string type =
"standard";
162 for (
int d=0; d < dim; d++) {
163 point.at(d) = vecPoints->at(i).at(d);
165 for (
int d=0; d < dofs; d++) {
166 result.at(d) = vecPoints->at(i).at(d);
168 vecBC_func_.at(loc)(&point.at(0), &result.at(0), t, &(vecBC_Parameters_.at(loc).at(0)));
170 for (UN j=0; j<blockMV->getNumVectors(); j++) {
171 Teuchos::ArrayRCP<SC> values = blockMV->getBlock(block)->getDataNonConst(j);
172 if ( !vecBCType_.at(loc).compare(
"Dirichlet") ) {
173 for (
int dd=0; dd < dofs; dd++)
174 values[ dofs * i + dd ] = result.at(dd);
176 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X") ) {
177 values[ dofs * i + 0 ] = result.at(0);
179 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y") ) {
180 values[ dofs * i + 1 ] = result.at(1);
182 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Z") ) {
183 values[ dofs * i + 2 ] = result.at(2);
185 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ) {
186 values[ dofs * i + 0 ] = result.at(0);
187 values[ dofs * i + 1 ] = result.at(1);
189 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ) {
190 values[ dofs * i + 0 ] = result.at(0);
191 values[ dofs * i + 2 ] = result.at(2);
193 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z") ) {
194 values[ dofs * i + 1 ] = result.at(1);
195 values[ dofs * i + 2 ] = result.at(2);
204 for (
int i=0; i<vecBCType_.size(); i++) {
205 if( vecBCType_.at(i) ==
"Neumann" && vecBlockID_.at(i)==block){
206 DomainPtr_Type domain = vecDomain_.at(i);
207 FEFacPtr_Type feFactory = Teuchos::rcp(
new FEFac_Type() );
208 feFactory->addFE(domain);
210 std::vector<SC> funcParameter(3 , 0.);
211 funcParameter[1] = t;
212 funcParameter[2] = vecFlag_.at(i);
213 int dim = domain->getDimension();
214 int dofs = vecDofs_.at(i);
215 MultiVectorPtr_Type a;
216 MultiVectorPtr_Type aUnique;
218 a = Teuchos::rcp(
new MultiVector_Type ( domain->getMapVecFieldRepeated() ) );
219 aUnique = Teuchos::rcp(
new MultiVector_Type ( domain->getMapVecFieldUnique() ) );
220 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a,
"Vector", vecBC_func_.at(i), funcParameter );
223 a = Teuchos::rcp(
new MultiVector_Type ( domain->getMapRepeated() ) );
224 aUnique = Teuchos::rcp(
new MultiVector_Type ( domain->getMapUnique() ) );
225 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a,
"Scalar", vecBC_func_.at(i), funcParameter );
227 aUnique->exportFromVector( a,
false,
"Add" );
228 MultiVectorPtr_Type mv = blockMV->getBlockNonConst(block);
229 mv->update(1.,*aUnique,1.);
235template<
class SC,
class LO,
class GO,
class NO>
238 Teuchos::ArrayRCP<SC> valuesEx = vecExternalSol_[loc]->getDataNonConst(0);
240 (*pointPtr_)[0] = valuesEx[index];
242 vecBC_func_.at(loc)( &((*pointPtr_)[0]), &((*resultPtr_)[0]), time, &(vecBC_Parameters_.at(loc).at(0)));
244 int dofs = resultPtr_->size();
245 if ( !vecBCType_.at(loc).compare(
"Dirichlet") ) {
246 if (type ==
"standard") {
247 for (
int dd=0; dd < dofs; dd++){
248 values[ dofs * index + dd ] = (*resultPtr_)[dd];
251 else if (type ==
"BCMinusVec"){
252 for (
int dd=0; dd < dofs; dd++)
253 values[ dofs * index + dd ] = (*resultPtr_)[dd] - valuesSubstract[ dofs * index + dd ];
255 else if (type ==
"VecMinusBC"){
256 for (
int dd=0; dd < dofs; dd++)
257 values[ dofs * index + dd ] = valuesSubstract[ dofs * index + dd ] - (*resultPtr_)[dd];
261 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
"BCBuilder::setDirichletBoundaryFromExternal only implemented for full Dirichlet.");
265template<
class SC,
class LO,
class GO,
class NO>
268 TEUCHOS_TEST_FOR_EXCEPTION( outBlockMV->getNumVectors()>1, std::runtime_error,
"BCBuilder::setRHS() only for getNumVectors == 1.");
274 for (
int block = 0; block < outBlockMV->size(); block++) {
276 if (vecFlowRateBool_[loc]) {
277 this->determineVelocityForFlowrate(loc,t);
279 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
280 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
281 int dim = vecDomain_.at(loc)->getDimension();
282 int dofs = vecDofs_.at(loc);
284 result.resize(dofs,0.);
285 point.resize(dim,0.);
287 resultPtr_ = Teuchos::rcp(
new vec_dbl_Type(dofs,0.));
288 pointPtr_ = Teuchos::rcp(
new vec_dbl_Type(dim,0.));
290 Teuchos::ArrayRCP<SC> valuesRHS = outBlockMV->getBlock(block)->getDataNonConst(0);
291 Teuchos::ArrayRCP<SC> valuesSubstract = substractBlockMV->getBlock(block)->getDataNonConst(0);
293 for (
int i = 0 ; i < bcFlags->size(); i++) {
294 int flag = bcFlags->at(i);
296 if (!vecBCType_.at(loc).compare(
"Dirichlet") ||
297 !vecBCType_.at(loc).compare(
"Dirichlet_X") ||
298 !vecBCType_.at(loc).compare(
"Dirichlet_Y") ||
299 !vecBCType_.at(loc).compare(
"Dirichlet_Z")||
300 !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ||
301 !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ||
302 !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z")) {
303 if (!vecExternalSol_[loc].is_null()) {
304 std::string type =
"BCMinusVec";
308 for (
int d=0; d < dim; d++) {
309 point.at(d) = vecPoints->at(i).at(d);
311 for (
int d=0; d < dofs; d++) {
312 result.at(d) = vecPoints->at(i).at(d);
314 vecBC_func_.at(loc)(&point.at(0), &result.at(0), t, &(vecBC_Parameters_.at(loc).at(0)));
316 for (UN j=0; j<outBlockMV->getNumVectors(); j++) {
317 Teuchos::ArrayRCP<SC> values = outBlockMV->getBlock(block)->getDataNonConst(j);
318 Teuchos::ArrayRCP<const SC> substractValues = substractBlockMV->getBlock(block)->getData(j);
319 if ( !vecBCType_.at(loc).compare(
"Dirichlet") ) {
320 for (
int dd=0; dd < dofs; dd++)
321 values[ dofs * i + dd ] = result.at(dd) - substractValues[ dofs * i + dd ];
323 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X") ) {
324 values[ dofs * i + 0 ] = result.at(0) - substractValues[ dofs * i + 0 ];
326 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y") ) {
327 values[ dofs * i + 1 ] = result.at(1) - substractValues[ dofs * i + 1 ];
329 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Z") ) {
330 values[ dofs * i + 2 ] = result.at(2) - substractValues[ dofs * i + 2 ];
332 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ) {
333 values[ dofs * i + 0 ] = result.at(0) - substractValues[ dofs * i + 0 ];
334 values[ dofs * i + 1 ] = result.at(1) - substractValues[ dofs * i + 1 ];
336 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ) {
337 values[ dofs * i + 0 ] = result.at(0) - substractValues[ dofs * i + 0 ];
338 values[ dofs * i + 2 ] = result.at(2) - substractValues[ dofs * i + 2 ];
340 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z") ) {
341 values[ dofs * i + 1 ] = result.at(1) - substractValues[ dofs * i + 1 ];
342 values[ dofs * i + 2 ] = result.at(2) - substractValues[ dofs * i + 2 ];
351 for (
int i=0; i<vecBCType_.size(); i++) {
352 if( vecBCType_.at(i) ==
"Neumann" && vecBlockID_.at(i)==block){
353 DomainPtr_Type domain = vecDomain_.at(i);
354 FEFacPtr_Type feFactory = Teuchos::rcp(
new FEFac_Type() );
355 feFactory->addFE(domain);
357 std::vector<SC> funcParameter(3 , 0.);
358 funcParameter[1] = t;
359 funcParameter[2] = vecFlag_.at(i);
360 int dim = domain->getDimension();
361 int dofs = vecDofs_.at(i);
362 MultiVectorPtr_Type a;
363 MultiVectorPtr_Type aUnique;
365 a = Teuchos::rcp(
new MultiVector_Type ( domain->getMapVecFieldRepeated() ) );
366 aUnique = Teuchos::rcp(
new MultiVector_Type ( domain->getMapVecFieldUnique() ) );
367 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a,
"Vector", vecBC_func_.at(i), funcParameter );
370 a = Teuchos::rcp(
new MultiVector_Type ( domain->getMapRepeated() ) );
371 aUnique = Teuchos::rcp(
new MultiVector_Type ( domain->getMapUnique() ) );
372 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a,
"Scalar", vecBC_func_.at(i), funcParameter );
374 aUnique->exportFromVector( a,
false,
"Add" );
375 MultiVectorPtr_Type mv = outBlockMV->getBlockNonConst(block);
376 mv->update(1.,*aUnique,1.);
377 MultiVectorPtr_Type mvMinus = substractBlockMV->getBlockNonConst(block);
378 mv->update(-1.,*mvMinus,1.);
384template<
class SC,
class LO,
class GO,
class NO>
390 TEUCHOS_TEST_FOR_EXCEPTION( outBlockMV->getNumVectors()>1, std::runtime_error,
"BCBuilder::setRHS() only for getNumVectors == 1.");
392 for (
int block = 0; block < outBlockMV->size(); block++) {
394 if (vecFlowRateBool_[loc]) {
395 this->determineVelocityForFlowrate(loc,t);
397 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
398 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
399 int dim = vecDomain_.at(loc)->getDimension();
400 int dofs = vecDofs_.at(loc);
402 result.resize(dofs,0.);
403 point.resize(dim,0.);
405 resultPtr_ = Teuchos::rcp(
new vec_dbl_Type(dofs,0.));
406 pointPtr_ = Teuchos::rcp(
new vec_dbl_Type(dim,0.));
408 Teuchos::ArrayRCP<SC> valuesRHS = outBlockMV->getBlock(block)->getDataNonConst(0);
409 Teuchos::ArrayRCP<SC> valuesSubstract = substractBlockMV->getBlock(block)->getDataNonConst(0);
410 for (
int i = 0 ; i < bcFlags->size(); i++) {
411 int flag = bcFlags->at(i);
413 if (!vecBCType_.at(loc).compare(
"Dirichlet") ||
414 !vecBCType_.at(loc).compare(
"Dirichlet_X") ||
415 !vecBCType_.at(loc).compare(
"Dirichlet_Y") ||
416 !vecBCType_.at(loc).compare(
"Dirichlet_Z") ||
417 !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ||
418 !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ||
419 !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z")) {
420 if (!vecExternalSol_[loc].is_null()) {
421 std::string type =
"VecMinusBC";
425 for (
int d=0; d < dim; d++) {
426 point.at(d) = vecPoints->at(i).at(d);
428 for (
int d=0; d < dofs; d++) {
429 result.at(d) = vecPoints->at(i).at(d);
431 vecBC_func_.at(loc)(&point.at(0), &result.at(0), t, &(vecBC_Parameters_.at(loc).at(0)));
433 for (UN j=0; j<outBlockMV->getNumVectors(); j++) {
434 Teuchos::ArrayRCP<SC> values = outBlockMV->getBlock(block)->getDataNonConst(j);
435 Teuchos::ArrayRCP<const SC> substractValues = substractBlockMV->getBlock(block)->getData(j);
436 if ( !vecBCType_.at(loc).compare(
"Dirichlet") ) {
437 for (
int dd=0; dd < dofs; dd++)
438 values[ dofs * i + dd ] = substractValues[ dofs * i + dd ] - result.at(dd);
440 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X") ) {
441 values[ dofs * i + 0 ] = substractValues[ dofs * i + 0 ] - result.at(0);
443 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y") ) {
444 values[ dofs * i + 1 ] = substractValues[ dofs * i + 1 ] - result.at(1);
446 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Z") ) {
447 values[ dofs * i + 2 ] = substractValues[ dofs * i + 2 ] - result.at(2);
449 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ) {
450 values[ dofs * i + 0 ] = substractValues[ dofs * i + 0 ] - result.at(0);
451 values[ dofs * i + 1 ] = substractValues[ dofs * i + 1 ] - result.at(1);
453 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ) {
454 values[ dofs * i + 0 ] = substractValues[ dofs * i + 0 ] - result.at(0);
455 values[ dofs * i + 2 ] = substractValues[ dofs * i + 2 ] - result.at(2);
457 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z") ) {
458 values[ dofs * i + 1 ] = substractValues[ dofs * i + 1 ] - result.at(1);
459 values[ dofs * i + 2 ] = substractValues[ dofs * i + 2 ] - result.at(2);
468 for (
int i=0; i<vecBCType_.size(); i++) {
469 if( vecBCType_.at(i) ==
"Neumann" && vecBlockID_.at(i)==block){
470 DomainPtr_Type domain = vecDomain_.at(i);
471 FEFacPtr_Type feFactory = Teuchos::rcp(
new FEFac_Type() );
472 feFactory->addFE(domain);
474 std::vector<SC> funcParameter(3 , 0.);
475 funcParameter[1] = t;
476 funcParameter[2] = vecFlag_.at(i);
477 int dim = domain->getDimension();
478 int dofs = vecDofs_.at(i);
479 MultiVectorPtr_Type a;
480 MultiVectorPtr_Type aUnique;
482 a = Teuchos::rcp(
new MultiVector_Type ( domain->getMapVecFieldRepeated() ) );
483 aUnique = Teuchos::rcp(
new MultiVector_Type ( domain->getMapVecFieldUnique() ) );
484 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a,
"Vector", vecBC_func_.at(i), funcParameter );
487 a = Teuchos::rcp(
new MultiVector_Type ( domain->getMapRepeated() ) );
488 aUnique = Teuchos::rcp(
new MultiVector_Type ( domain->getMapUnique() ) );
489 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a,
"Scalar", vecBC_func_.at(i), funcParameter );
491 aUnique->exportFromVector( a,
false,
"Add" );
492 MultiVectorPtr_Type mv = outBlockMV->getBlockNonConst(block);
493 mv->update(1.,*aUnique,1.);
494 MultiVectorPtr_Type mvMinus = substractBlockMV->getBlockNonConst(block);
495 mv->update(1.,*mvMinus,-1.);
501template<
class SC,
class LO,
class GO,
class NO>
505 for (
int block = 0; block < blockMV->size(); block++) {
507 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
508 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
509 int dim = vecDomain_.at(loc)->getDimension();
510 int dofs = vecDofs_.at(loc);
512 result.resize(dofs,0.);
514 for (
int i = 0 ; i < bcFlags->size(); i++) {
515 int flag = bcFlags->at(i);
517 if ( !vecBCType_.at(loc).compare(
"Dirichlet") ||
518 !vecBCType_.at(loc).compare(
"Dirichlet_X") ||
519 !vecBCType_.at(loc).compare(
"Dirichlet_Y") ||
520 !vecBCType_.at(loc).compare(
"Dirichlet_Z") ||
521 !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ||
522 !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ||
523 !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z")
526 for (UN j=0; j<blockMV->getNumVectors(); j++) {
527 Teuchos::ArrayRCP<SC> values = blockMV->getBlock(block)->getDataNonConst(j);
528 if ( !vecBCType_.at(loc).compare(
"Dirichlet") ) {
529 for (
int dd=0; dd < dofs; dd++)
530 values[ dofs * i + dd ] = result[dd];
532 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X") ) {
533 values[ dofs * i + 0 ] = result[0];
535 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y") ) {
536 values[ dofs * i + 1 ] = result[1];
538 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Z") ) {
539 values[ dofs * i + 2 ] = result[2];
541 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ) {
542 values[ dofs * i + 0 ] = result[0];
543 values[ dofs * i + 1 ] = result[1];
545 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ) {
546 values[ dofs * i + 0 ] = result[0];
547 values[ dofs * i + 2 ] = result[2];
549 else if ( !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z") ) {
550 values[ dofs * i + 1 ] = result[1];
551 values[ dofs * i + 2 ] = result[2];
572template<
class SC,
class LO,
class GO,
class NO>
573void BCBuilder<SC,LO,GO,NO>::determineVelocityForFlowrate(LO i,
double time)
const{
576 DomainPtr_Type domain = vecDomain_.at(i);
578 MultiVectorConstPtr_Type parabolic_unique_const = vecExternalSol_[i];
579 MultiVectorPtr_Type parabolic_unique = Teuchos::rcp_const_cast<MultiVector_Type> ( parabolic_unique_const );
581 MultiVectorPtr_Type parabolic_rep = Teuchos::rcp(
new MultiVector_Type ( domain->getMapRepeated() ) );
583 SC maxValue = parabolic_unique->getMax();
584 parabolic_unique->scale(1./maxValue);
585 parabolic_rep->importFromVector(parabolic_unique,
false,
"Insert");
587 FEFacPtr_Type feFactory = Teuchos::rcp(
new FEFac_Type() );
588 feFactory->addFE(domain);
590 std::vector<SC> funcParameter = vecBC_Parameters_[i];
591 vec_dbl_Type p1 = {1.,1.,1.};
592 vec_dbl_Type flowRate = {0.};
593 vecBC_func_flowRate_.at(i)( &(p1[0]), &(flowRate[0]), time, &(funcParameter[0]));
595 double flowRateParabolic=0.;
596 feFactory->assemblyFlowRate(domain->getDimension(), flowRateParabolic, domain->getFEType(),1, vecFlag_[i] , parabolic_rep);
598 double maxVelocity = flowRate[0] / std::fabs(flowRateParabolic);
600 vecBC_Parameters_[i][0] = maxVelocity;
606template<
class SC,
class LO,
class GO,
class NO>
613template<
class SC,
class LO,
class GO,
class NO>
619 std::vector<int>::const_iterator it = vecBlockID_.begin();
620 while (it!=vecBlockID_.end() && !found) {
621 it = std::find(it,vecBlockID_.end(),block);
622 if (it!=vecBlockID_.end()) {
623 loc = distance(vecBlockID_.begin(),it);
625 if (!vecBCType_.at(loc).compare(
"Dirichlet") ||
626 !vecBCType_.at(loc).compare(
"Dirichlet_X") ||
627 !vecBCType_.at(loc).compare(
"Dirichlet_Y") ||
628 !vecBCType_.at(loc).compare(
"Dirichlet_Z") ||
629 !vecBCType_.at(loc).compare(
"Dirichlet_X_Y") ||
630 !vecBCType_.at(loc).compare(
"Dirichlet_X_Z") ||
631 !vecBCType_.at(loc).compare(
"Dirichlet_Y_Z")
642template<
class SC,
class LO,
class GO,
class NO>
645#ifdef BCBuilder_TIMER
646 Teuchos::TimeMonitor FindFlagMonitor(*FindFlagTimer_);
649 bool hasFlag =
false;
652 std::vector<int>::const_iterator it = vecFlag_.begin();
653 while (it!=vecFlag_.end() && !found) {
654 it = std::find(it,vecFlag_.end(),flag);
655 if (it!=vecFlag_.end()) {
656 loc = distance(vecFlag_.begin(),it);
657 if (vecBlockID_.at(loc)==block) {
671template<
class SC,
class LO,
class GO,
class NO>
674 UN numBlocks = blockMatrix->size();
676 for (UN blockRow = 0; blockRow < numBlocks; blockRow++) {
678#ifdef BCBuilder_TIMER
679 Teuchos::TimeMonitor SetSystemRowMonitor(*SetSystemRowTimer_);
681 bool boolBlockHasDirichlet =
false;
684#ifdef BCBuilder_TIMER
685 Teuchos::TimeMonitor BlockRowHasDirichletMonitor(*BlockRowHasDirichletTimer_);
690 if (boolBlockHasDirichlet){
691 for (UN blockCol = 0; blockCol < numBlocks ; blockCol++) {
692 if ( blockMatrix->blockExists( blockRow, blockCol ) ) {
693 MatrixPtr_Type matrix = blockMatrix->getBlock( blockRow, blockCol );
702template<
class SC,
class LO,
class GO,
class NO>
705 UN numBlocks = blockMatrix->size();
707 for (UN blockRow = 0; blockRow < numBlocks; blockRow++) {
709#ifdef BCBuilder_TIMER
710 Teuchos::TimeMonitor SetSystemRowMonitor(*SetSystemRowTimer_);
712 bool boolBlockHasDirichlet =
false;
715#ifdef BCBuilder_TIMER
716 Teuchos::TimeMonitor BlockRowHasDirichletMonitor(*BlockRowHasDirichletTimer_);
722 if (boolBlockHasDirichlet){
723 for (UN blockCol = 0; blockCol < numBlocks ; blockCol++) {
724 if (blockMatrix->blockExists(blockRow, blockCol)) {
725 MatrixPtr_Type matrix = blockMatrix->getBlock(blockRow, blockCol);
727 }
else if (blockRow == blockCol) {
734 auto matrixMap = this->vecDomain_.at(loc)->getMapUnique();
736 auto tpetraMatrix = Teuchos::rcp(
new Tpetra::CrsMatrix<SC,LO,GO,NO>(matrixMap->getTpetraMap(), matrixMap->getTpetraMap(), 1));
738 MatrixPtr_Type matrix = Teuchos::rcp(
new Matrix_Type(tpetraMatrix));
741 Teuchos::Array<LO> colIndex(1, 0);
742 Teuchos::Array<SC> val(1, Teuchos::ScalarTraits<SC>::zero());
743 for (
auto i = 0; i < matrixMap->getNodeNumElements(); i++) {
746 matrix->insertLocalValues(i, colIndex, val);
748 matrix->fillComplete(matrixMap, matrixMap);
750 blockMatrix->addBlock(matrix, blockRow, blockCol);
758template<
class SC,
class LO,
class GO,
class NO>
761 matrix->resumeFill();
763 UN dofsPerNode = vecDofs_.at(loc);
764 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
765 MapConstPtr_Type nodeMap = vecDomain_.at(loc)->getMapUnique();
767 for (LO i=0; i<bcFlags->size(); i++) {
768 int flag = bcFlags->at(i);
772 if(
findFlag(flag, blockRow, locThisFlag) ){
774 if( !vecBCType_.at(locThisFlag).compare(
"Dirichlet") ||
775 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_X") ||
776 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_Y") ||
777 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_Z") ||
778 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_X_Y") ||
779 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_X_Z") ||
780 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_Y_Z") )
791 matrix->fillComplete( matrix->getMap(
"domain"), matrix->getMap(
"range") );
795template<
class SC,
class LO,
class GO,
class NO>
798 Teuchos::ArrayView<const SC> valuesOld;
799 Teuchos::ArrayView<const LO> indices;
801 MapConstPtr_Type colMap = matrix->getMap(
"col");
802 LO localDof = (LO) (dofsPerNode * localNode);
803 for (UN dof=0; dof<dofsPerNode; dof++) {
804 if ( vecBCType_.at(loc) ==
"Dirichlet" ||
805 (vecBCType_.at(loc) ==
"Dirichlet_X" && dof==0 ) ||
806 (vecBCType_.at(loc) ==
"Dirichlet_Y" && dof==1 ) ||
807 (vecBCType_.at(loc) ==
"Dirichlet_Z" && dof==2 ) ||
808 (vecBCType_.at(loc) ==
"Dirichlet_X_Y" && dof!=2 )||
809 (vecBCType_.at(loc) ==
"Dirichlet_X_Z" && dof!=1 )||
810 (vecBCType_.at(loc) ==
"Dirichlet_Y_Z" && dof!=0 )
813 GO globalDof = matrix->getMap()->getGlobalElement( localDof );
814 matrix->getLocalRowView(localDof, indices, valuesOld);
816 for (UN j=0; j<indices.size(); j++) {
817 rowSum += abs(valuesOld[j]);
819 Teuchos::Array<SC> values( valuesOld.size(), Teuchos::ScalarTraits<SC>::zero() );
821 for (UN j=0; j<indices.size() && !setOne; j++) {
822 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
823 values[j] = valuesOld[j]*eps;
828 matrix->replaceLocalValues(localDof, indices(), values());
836template<
class SC,
class LO,
class GO,
class NO>
839 matrix->resumeFill();
841 UN dofsPerNode = vecDofs_.at(loc);
842 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
844 for (LO i=0; i<bcFlags->size(); i++) {
845 int flag = bcFlags->at(i);
849 if(
findFlag(flag, blockRow, locThisFlag) ){
851 if( !vecBCType_.at(locThisFlag).compare(
"Dirichlet") ||
852 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_X") ||
853 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_Y") ||
854 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_Z") ||
855 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_X_Y") ||
856 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_X_Z") ||
857 !vecBCType_.at(locThisFlag).compare(
"Dirichlet_Y_Z") )
868 matrix->fillComplete( matrix->getMap(
"domain"), matrix->getMap(
"range") );
871template<
class SC,
class LO,
class GO,
class NO>
874 Teuchos::ArrayView<const SC> valuesOld;
875 Teuchos::ArrayView<const LO> indices;
877 MapConstPtr_Type colMap = matrix->getMap(
"col");
878 LO localDof = (LO) (dofsPerNode * localNode);
879 for (UN dof=0; dof<dofsPerNode; dof++) {
880 if ( vecBCType_.at(loc) ==
"Dirichlet" ||
881 (vecBCType_.at(loc) ==
"Dirichlet_X" && dof==0 ) ||
882 (vecBCType_.at(loc) ==
"Dirichlet_Y" && dof==1 ) ||
883 (vecBCType_.at(loc) ==
"Dirichlet_Z" && dof==2 ) ||
884 (vecBCType_.at(loc) ==
"Dirichlet_X_Y" && dof!=2 )||
885 (vecBCType_.at(loc) ==
"Dirichlet_X_Z" && dof!=1 )||
886 (vecBCType_.at(loc) ==
"Dirichlet_Y_Z" && dof!=0 )
889 GO globalDof = matrix->getMap()->getGlobalElement( localDof );
890 matrix->getLocalRowView(localDof, indices, valuesOld);
891 Teuchos::Array<SC> values( valuesOld.size(), Teuchos::ScalarTraits<SC>::zero() );
893 for (UN j=0; j<indices.size() && !setOne; j++) {
894 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
895 values[j] = Teuchos::ScalarTraits<SC>::one();
899 matrix->replaceLocalValues(localDof, indices(), values());
906template<
class SC,
class LO,
class GO,
class NO>
909 Teuchos::ArrayView<const SC> valuesOld;
910 Teuchos::ArrayView<const LO> indices;
911 LO localDof = (LO) (dofsPerNode * localNode);
912 for (UN dof=0; dof<dofsPerNode; dof++) {
913 if ( vecBCType_.at(loc) ==
"Dirichlet" ||
914 (vecBCType_.at(loc) ==
"Dirichlet_X" && dof==0 ) ||
915 (vecBCType_.at(loc) ==
"Dirichlet_Y" && dof==1 ) ||
916 (vecBCType_.at(loc) ==
"Dirichlet_Z" && dof==2 ) ||
917 (vecBCType_.at(loc) ==
"Dirichlet_X_Y" && dof!=2 )||
918 (vecBCType_.at(loc) ==
"Dirichlet_X_Z" && dof!=1 )||
919 (vecBCType_.at(loc) ==
"Dirichlet_Y_Z" && dof!=0 )
921 matrix->getLocalRowView(localDof, indices, valuesOld);
922 Teuchos::Array<SC> values( valuesOld.size(), Teuchos::ScalarTraits<SC>::zero() );
923 matrix->replaceLocalValues(localDof, indices(), values());
931template<
class SC,
class LO,
class GO,
class NO>
933 std::vector<int>::iterator it;
934 it = find(vecBlockID_.begin(),vecBlockID_.end(),block);
935#ifdef ASSERTS_WARNINGS
936 MYASSERT(it!=vecBlockID_.end(),
"Requested dof information is not known to BCBuilder class.");
939 return vecDofs_.at(distance(vecBlockID_.begin(),it));
void setBCMinusVector(const BlockMultiVectorPtr_Type &outBlockMV, const BlockMultiVectorPtr_Type &substractBlockMV, double t=0.) const
Definition BCBuilder_def.hpp:266
void setVectorMinusBC(const BlockMultiVectorPtr_Type &outBlockMV, const BlockMultiVectorPtr_Type &substractBlockMV, double t=0.) const
Definition BCBuilder_def.hpp:385
void setRHS(const BlockMultiVectorPtr_Type &blockMV, double t=0.) const
Setting boundary conditions to (block)vector.
Definition BCBuilder_def.hpp:123
void set(const BlockMatrixPtr_Type &blockMatrix, const BlockMultiVectorPtr_Type &b, double t=0.) const
Setting bundary condtions to problem.
Definition BCBuilder_def.hpp:114
void setSystem(const BlockMatrixPtr_Type &blockMatrix) const
Set boundary conditions to system.
Definition BCBuilder_def.hpp:703
void addBC(BC_func_Type funcBC, int flag, int block, const DomainPtr_Type &domain, std::string type, int dofs)
Adding Boundary Condition.
Definition BCBuilder_def.hpp:46
void setLocalRowZero(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc) const
Definition BCBuilder_def.hpp:907
void setDirichletBoundaryFromExternal(Teuchos::ArrayRCP< SC > &values, LO index, int loc, double time, std::string type, Teuchos::ArrayRCP< SC > valuesSubstract=Teuchos::null) const
Definition BCBuilder_def.hpp:236
void setSystemScaled(const BlockMatrixPtr_Type &blockMatrix, double eps=1.0) const
Set boundary conditions to system.
Definition BCBuilder_def.hpp:672
void setDirichletBC(const MatrixPtr_Type &matrix, int loc, int blockRow, bool isDiagonalBlock) const
Definition BCBuilder_def.hpp:837
void setAllDirichletZero(const BlockMultiVectorPtr_Type &blockMV) const
Definition BCBuilder_def.hpp:502
bool blockHasDirichletBC(int block) const
Definition BCBuilder_def.hpp:607
void setDirichletBCScaled(const MatrixPtr_Type &matrix, int loc, int blockRow, bool isDiagonalBlock, double eps=1.0) const
Definition BCBuilder_def.hpp:759
bool findFlag(LO flag, int block, int &loc) const
Definition BCBuilder_def.hpp:643
void setLocalRowEntry(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc, double eps) const
Definition BCBuilder_def.hpp:796
void setLocalRowOne(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc) const
Definition BCBuilder_def.hpp:872
int dofsPerNodeAtBlock(int block)
Definition BCBuilder_def.hpp:932
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36