Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
BCBuilder_def.hpp
1#ifndef BCBuilder_def_hpp
2#define BCBuilder_def_hpp
3
4#include "BCBuilder_decl.hpp"
5#include "feddlib/core/Utils/FEDDUtils.hpp"
6#include <Teuchos_RCPDecl.hpp>
7#include <Teuchos_ScalarTraitsDecl.hpp>
8#include <Xpetra_MatrixFactory.hpp>
9
18
19namespace FEDD {
20template<class SC,class LO,class GO,class NO>
21BCBuilder<SC,LO,GO,NO>::BCBuilder():
22 vecBC_func_(),
23 vecFlag_(),
24 vecBlockID_(),
25 vecDomain_(),
26 vecBCType_(),
27 vecDofs_(),
28 vecBC_Parameters_(),
29 vecExternalSol_(0),
30 resultPtr_(),
31 pointPtr_()
32#ifdef BCBuilder_TIMER
33,SetSystemRowTimer_(Teuchos::TimeMonitor::getNewCounter("BCBuilder: SetSystemRow")),
34BlockRowHasDirichletTimer_(Teuchos::TimeMonitor::getNewCounter("BCBuilder: SetSystemRow: BlockHasDiri")),
35FindFlagTimer_(Teuchos::TimeMonitor::getNewCounter("BCBuilder: SetSystemRow: BlockHasDiri: FindFlag"))
36#endif
37{
38
39}
40
41template<class SC,class LO,class GO,class NO>
42void BCBuilder<SC,LO,GO,NO>::addBC(BC_func_Type funcBC, int flag, int block, const DomainPtr_Type &domain, std::string type, int dofs){
43
44 vecBC_func_.push_back(funcBC);
45 vecFlag_.push_back(flag);
46 vecBlockID_.push_back(block);
47 vecDomain_.push_back(domain);
48 vecBCType_.push_back(type);
49 vecDofs_.push_back(dofs);
50 vec_dbl_Type dummy(1,0.);
51 vecBC_Parameters_.push_back(dummy);
52 MultiVectorConstPtr_Type dummyMV;
53 vecExternalSol_.push_back(dummyMV);
54
55}
56
57template<class SC,class LO,class GO,class NO>
58void 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 &parameter_vec){
59
60 vecBC_func_.push_back(funcBC);
61 vecFlag_.push_back(flag);
62 vecBlockID_.push_back(block);
63 vecDomain_.push_back(domain);
64 vecBCType_.push_back(type);
65 vecDofs_.push_back(dofs);
66 vecBC_Parameters_.push_back(parameter_vec);
67 MultiVectorConstPtr_Type dummyMV;
68 vecExternalSol_.push_back(dummyMV);
69
70}
71
72template<class SC,class LO,class GO,class NO>
73void 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 &parameter_vec, MultiVectorConstPtr_Type& externalSol){
74
75 vecBC_func_.push_back(funcBC);
76 vecFlag_.push_back(flag);
77 vecBlockID_.push_back(block);
78 vecDomain_.push_back(domain);
79 vecBCType_.push_back(type);
80 vecDofs_.push_back(dofs);
81 vecBC_Parameters_.push_back(parameter_vec);
82 vecExternalSol_.push_back(externalSol);
83
84}
85
86
87template<class SC,class LO,class GO,class NO>
88void BCBuilder<SC,LO,GO,NO>::set(const BlockMatrixPtr_Type &blockMatrix, const BlockMultiVectorPtr_Type &blockMV, double t) const{
89
90 setSystem(blockMatrix);
91
92 setRHS(blockMV, t);
93
94}
95
96template<class SC,class LO,class GO,class NO>
97void BCBuilder<SC,LO,GO,NO>::setRHS(const BlockMultiVectorPtr_Type &blockMV, double t) const{
98 vec_dbl_Type result;
99 vec_dbl_Type point;
100
101 TEUCHOS_TEST_FOR_EXCEPTION( blockMV->getNumVectors()>1, std::runtime_error, "BCBuilder::setRHS() only for getNumVectors == 1.");
102
103 int loc = 0;
104 for (int block = 0; block < blockMV->size(); block++) { // blocks of RHS vector
105 if(blockHasDirichletBC(block, loc)){
106
107 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
108 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
109 int dim = vecDomain_.at(loc)->getDimension();
110 int dofs = vecDofs_.at(loc);
111
112 result.resize(dofs,0.);
113 point.resize(dim,0.);
114
115 resultPtr_ = Teuchos::rcp(new vec_dbl_Type(dofs,0.));
116 pointPtr_ = Teuchos::rcp(new vec_dbl_Type(dim,0.));
117
118 Teuchos::ArrayRCP<SC> valuesRHS = blockMV->getBlock(block)->getDataNonConst(0);
119 for (int i = 0 ; i < bcFlags->size(); i++) {
120 int flag = bcFlags->at(i);
121 if(findFlag(flag, block, loc)){ // loc changed to postion corresponding to bcFlag
122 if (!vecBCType_.at(loc).compare("Dirichlet") ||
123 !vecBCType_.at(loc).compare("Dirichlet_X") ||
124 !vecBCType_.at(loc).compare("Dirichlet_Y") ||
125 !vecBCType_.at(loc).compare("Dirichlet_Z")||
126 !vecBCType_.at(loc).compare("Dirichlet_X_Y") ||
127 !vecBCType_.at(loc).compare("Dirichlet_X_Z") ||
128 !vecBCType_.at(loc).compare("Dirichlet_Y_Z")) {
129 if (!vecExternalSol_[loc].is_null()) { // we use the external vector here
130 std::string type = "standard";
131 this->setDirichletBoundaryFromExternal(valuesRHS, i, loc, t, type);
132 }
133 else {
134 for (int d=0; d < dim; d++) {
135 point.at(d) = vecPoints->at(i).at(d);
136 }
137 for (int d=0; d < dofs; d++) {
138 result.at(d) = vecPoints->at(i).at(d);
139 }
140 vecBC_func_.at(loc)(&point.at(0), &result.at(0), t, &(vecBC_Parameters_.at(loc).at(0)));
141
142 for (UN j=0; j<blockMV->getNumVectors(); j++) {
143 Teuchos::ArrayRCP<SC> values = blockMV->getBlock(block)->getDataNonConst(j);
144 if ( !vecBCType_.at(loc).compare("Dirichlet") ) {
145 for (int dd=0; dd < dofs; dd++)
146 values[ dofs * i + dd ] = result.at(dd);
147 }
148 else if ( !vecBCType_.at(loc).compare("Dirichlet_X") ) {
149 values[ dofs * i + 0 ] = result.at(0);
150 }
151 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y") ) {
152 values[ dofs * i + 1 ] = result.at(1);
153 }
154 else if ( !vecBCType_.at(loc).compare("Dirichlet_Z") ) {
155 values[ dofs * i + 2 ] = result.at(2);
156 }
157 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Y") ) {
158 values[ dofs * i + 0 ] = result.at(0);
159 values[ dofs * i + 1 ] = result.at(1);
160 }
161 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Z") ) {
162 values[ dofs * i + 0 ] = result.at(0);
163 values[ dofs * i + 2 ] = result.at(2);
164 }
165 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y_Z") ) {
166 values[ dofs * i + 1 ] = result.at(1);
167 values[ dofs * i + 2 ] = result.at(2);
168 }
169 }
170 }
171 }
172 }
173 }
174 }
175 //Code below is experimental
176 for (int i=0; i<vecBCType_.size(); i++) {
177 if( vecBCType_.at(i) == "Neumann" && vecBlockID_.at(i)==block){
178 DomainPtr_Type domain = vecDomain_.at(i);
179 FEFacPtr_Type feFactory = Teuchos::rcp( new FEFac_Type() );
180 feFactory->addFE(domain);
181
182 std::vector<SC> funcParameter(3 , 0.);//0: order, 1:time, 2:surface flag (no specific flag is used in function above)
183 funcParameter[1] = t;
184 funcParameter[2] = vecFlag_.at(i);
185 int dim = domain->getDimension();
186 int dofs = vecDofs_.at(i);
187 MultiVectorPtr_Type a;
188 MultiVectorPtr_Type aUnique;
189 if (dofs>1){
190 a = Teuchos::rcp(new MultiVector_Type ( domain->getMapVecFieldRepeated() ) );
191 aUnique = Teuchos::rcp(new MultiVector_Type ( domain->getMapVecFieldUnique() ) );
192 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a, "Vector", vecBC_func_.at(i), funcParameter );
193 }
194 else{
195 a = Teuchos::rcp(new MultiVector_Type ( domain->getMapRepeated() ) );
196 aUnique = Teuchos::rcp(new MultiVector_Type ( domain->getMapUnique() ) );
197 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a, "Scalar", vecBC_func_.at(i), funcParameter );
198 }
199 aUnique->exportFromVector( a, false, "Add" );
200 MultiVectorPtr_Type mv = blockMV->getBlockNonConst(block);
201 mv->update(1.,*aUnique,1.);
202 }
203 }
204 }
205}
206
207template<class SC,class LO,class GO,class NO>
208void BCBuilder<SC,LO,GO,NO>::setDirichletBoundaryFromExternal(Teuchos::ArrayRCP<SC>& values/*values will be set to this vector*/, LO index, int loc, double time, std::string type, Teuchos::ArrayRCP<SC> valuesSubstract) const{
209
210 Teuchos::ArrayRCP<SC> valuesEx = vecExternalSol_[loc]->getDataNonConst(0);//We assume that our Multivector only has one column
211
212 (*pointPtr_)[0] = valuesEx[index]; // laplace solution in current point
213
214 vecBC_func_.at(loc)( &((*pointPtr_)[0]), &((*resultPtr_)[0]), time, &(vecBC_Parameters_.at(loc).at(0)));
215
216 int dofs = resultPtr_->size();
217 if ( !vecBCType_.at(loc).compare("Dirichlet") ) {
218 if (type == "standard") {
219 for (int dd=0; dd < dofs; dd++){
220 values[ dofs * index + dd ] = (*resultPtr_)[dd];
221 }
222 }
223 else if (type == "BCMinusVec"){
224 for (int dd=0; dd < dofs; dd++)
225 values[ dofs * index + dd ] = (*resultPtr_)[dd] - valuesSubstract[ dofs * index + dd ];
226 }
227 else if (type == "VecMinusBC"){
228 for (int dd=0; dd < dofs; dd++)
229 values[ dofs * index + dd ] = valuesSubstract[ dofs * index + dd ] - (*resultPtr_)[dd];
230 }
231 }
232 else {
233 TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, "BCBuilder::setDirichletBoundaryFromExternal only implemented for full Dirichlet.");
234 }
235}
236
237template<class SC,class LO,class GO,class NO>
238void BCBuilder<SC,LO,GO,NO>::setBCMinusVector(const BlockMultiVectorPtr_Type &outBlockMV, const BlockMultiVectorPtr_Type &substractBlockMV, double t ) const{
239
240 TEUCHOS_TEST_FOR_EXCEPTION( outBlockMV->getNumVectors()>1, std::runtime_error, "BCBuilder::setRHS() only for getNumVectors == 1.");
241
242 vec_dbl_Type result;
243 vec_dbl_Type point;
244 int loc = 0;
245
246 for (int block = 0; block < outBlockMV->size(); block++) { // blocks of RHS vector
247 if(blockHasDirichletBC(block, loc)){
248 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
249 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
250 int dim = vecDomain_.at(loc)->getDimension();
251 int dofs = vecDofs_.at(loc);
252
253 result.resize(dofs,0.);
254 point.resize(dim,0.);
255
256 resultPtr_ = Teuchos::rcp(new vec_dbl_Type(dofs,0.));
257 pointPtr_ = Teuchos::rcp(new vec_dbl_Type(dim,0.));
258
259 Teuchos::ArrayRCP<SC> valuesRHS = outBlockMV->getBlock(block)->getDataNonConst(0);
260 Teuchos::ArrayRCP<SC> valuesSubstract = substractBlockMV->getBlock(block)->getDataNonConst(0);
261
262 for (int i = 0 ; i < bcFlags->size(); i++) {
263 int flag = bcFlags->at(i);
264 if(findFlag(flag,block, loc)){ // loc changed to postion corresponding to bcFlag
265 if (!vecBCType_.at(loc).compare("Dirichlet") ||
266 !vecBCType_.at(loc).compare("Dirichlet_X") ||
267 !vecBCType_.at(loc).compare("Dirichlet_Y") ||
268 !vecBCType_.at(loc).compare("Dirichlet_Z")||
269 !vecBCType_.at(loc).compare("Dirichlet_X_Y") ||
270 !vecBCType_.at(loc).compare("Dirichlet_X_Z") ||
271 !vecBCType_.at(loc).compare("Dirichlet_Y_Z")) {
272 if (!vecExternalSol_[loc].is_null()) { // we use the external vector here
273 std::string type = "BCMinusVec";
274 this->setDirichletBoundaryFromExternal(valuesRHS, i, loc, t, type, valuesSubstract);
275 }
276 else{
277 for (int d=0; d < dim; d++) {
278 point.at(d) = vecPoints->at(i).at(d);
279 }
280 for (int d=0; d < dofs; d++) {
281 result.at(d) = vecPoints->at(i).at(d);
282 }
283 vecBC_func_.at(loc)(&point.at(0), &result.at(0), t, &(vecBC_Parameters_.at(loc).at(0)));
284
285 for (UN j=0; j<outBlockMV->getNumVectors(); j++) {
286 Teuchos::ArrayRCP<SC> values = outBlockMV->getBlock(block)->getDataNonConst(j);
287 Teuchos::ArrayRCP<const SC> substractValues = substractBlockMV->getBlock(block)->getData(j);
288 if ( !vecBCType_.at(loc).compare("Dirichlet") ) {
289 for (int dd=0; dd < dofs; dd++)
290 values[ dofs * i + dd ] = result.at(dd) - substractValues[ dofs * i + dd ];
291 }
292 else if ( !vecBCType_.at(loc).compare("Dirichlet_X") ) {
293 values[ dofs * i + 0 ] = result.at(0) - substractValues[ dofs * i + 0 ];
294 }
295 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y") ) {
296 values[ dofs * i + 1 ] = result.at(1) - substractValues[ dofs * i + 1 ];
297 }
298 else if ( !vecBCType_.at(loc).compare("Dirichlet_Z") ) {
299 values[ dofs * i + 2 ] = result.at(2) - substractValues[ dofs * i + 2 ];
300 }
301 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Y") ) {
302 values[ dofs * i + 0 ] = result.at(0) - substractValues[ dofs * i + 0 ];
303 values[ dofs * i + 1 ] = result.at(1) - substractValues[ dofs * i + 1 ];
304 }
305 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Z") ) {
306 values[ dofs * i + 0 ] = result.at(0) - substractValues[ dofs * i + 0 ];
307 values[ dofs * i + 2 ] = result.at(2) - substractValues[ dofs * i + 2 ];
308 }
309 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y_Z") ) {
310 values[ dofs * i + 1 ] = result.at(1) - substractValues[ dofs * i + 1 ];
311 values[ dofs * i + 2 ] = result.at(2) - substractValues[ dofs * i + 2 ];
312 }
313 }
314 }
315 }
316 }
317 }
318 }
319 //Code below is experimental
320 for (int i=0; i<vecBCType_.size(); i++) {
321 if( vecBCType_.at(i) == "Neumann" && vecBlockID_.at(i)==block){
322 DomainPtr_Type domain = vecDomain_.at(i);
323 FEFacPtr_Type feFactory = Teuchos::rcp( new FEFac_Type() );
324 feFactory->addFE(domain);
325
326 std::vector<SC> funcParameter(3 , 0.);//0: order, 1:time, 2:surface flag (no specific flag is used in function above)
327 funcParameter[1] = t;
328 funcParameter[2] = vecFlag_.at(i);
329 int dim = domain->getDimension();
330 int dofs = vecDofs_.at(i);
331 MultiVectorPtr_Type a;
332 MultiVectorPtr_Type aUnique;
333 if (dofs>1){
334 a = Teuchos::rcp(new MultiVector_Type ( domain->getMapVecFieldRepeated() ) );
335 aUnique = Teuchos::rcp(new MultiVector_Type ( domain->getMapVecFieldUnique() ) );
336 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a, "Vector", vecBC_func_.at(i), funcParameter );
337 }
338 else{
339 a = Teuchos::rcp(new MultiVector_Type ( domain->getMapRepeated() ) );
340 aUnique = Teuchos::rcp(new MultiVector_Type ( domain->getMapUnique() ) );
341 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a, "Scalar", vecBC_func_.at(i), funcParameter );
342 }
343 aUnique->exportFromVector( a, false, "Add" );
344 MultiVectorPtr_Type mv = outBlockMV->getBlockNonConst(block);
345 mv->update(1.,*aUnique,1.);
346 MultiVectorPtr_Type mvMinus = substractBlockMV->getBlockNonConst(block);
347 mv->update(-1.,*mvMinus,1.);
348 }
349 }
350 }
351}
352
353template<class SC,class LO,class GO,class NO>
354void BCBuilder<SC,LO,GO,NO>::setVectorMinusBC(const BlockMultiVectorPtr_Type &outBlockMV, const BlockMultiVectorPtr_Type &substractBlockMV, double t) const{
355 vec_dbl_Type result;
356 vec_dbl_Type point;
357 int loc = 0;
358
359 TEUCHOS_TEST_FOR_EXCEPTION( outBlockMV->getNumVectors()>1, std::runtime_error, "BCBuilder::setRHS() only for getNumVectors == 1.");
360
361 for (int block = 0; block < outBlockMV->size(); block++) { // blocks of RHS vector
362 if(blockHasDirichletBC(block, loc)){
363 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
364 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
365 int dim = vecDomain_.at(loc)->getDimension();
366 int dofs = vecDofs_.at(loc);
367
368 result.resize(dofs,0.);
369 point.resize(dim,0.);
370
371 resultPtr_ = Teuchos::rcp(new vec_dbl_Type(dofs,0.));
372 pointPtr_ = Teuchos::rcp(new vec_dbl_Type(dim,0.));
373
374 Teuchos::ArrayRCP<SC> valuesRHS = outBlockMV->getBlock(block)->getDataNonConst(0);
375 Teuchos::ArrayRCP<SC> valuesSubstract = substractBlockMV->getBlock(block)->getDataNonConst(0);
376 for (int i = 0 ; i < bcFlags->size(); i++) {
377 int flag = bcFlags->at(i);
378 if(findFlag(flag,block, loc)){ // loc changed to postion corresponding to bcFlag
379 if (!vecBCType_.at(loc).compare("Dirichlet") ||
380 !vecBCType_.at(loc).compare("Dirichlet_X") ||
381 !vecBCType_.at(loc).compare("Dirichlet_Y") ||
382 !vecBCType_.at(loc).compare("Dirichlet_Z") ||
383 !vecBCType_.at(loc).compare("Dirichlet_X_Y") ||
384 !vecBCType_.at(loc).compare("Dirichlet_X_Z") ||
385 !vecBCType_.at(loc).compare("Dirichlet_Y_Z")) {
386 if (!vecExternalSol_[loc].is_null()) { // we use the external vector here
387 std::string type = "VecMinusBC";
388 this->setDirichletBoundaryFromExternal(valuesRHS, i, loc, t, type, valuesSubstract);
389 }
390 else{
391 for (int d=0; d < dim; d++) {
392 point.at(d) = vecPoints->at(i).at(d);
393 }
394 for (int d=0; d < dofs; d++) {
395 result.at(d) = vecPoints->at(i).at(d);
396 }
397 vecBC_func_.at(loc)(&point.at(0), &result.at(0), t, &(vecBC_Parameters_.at(loc).at(0)));
398
399 for (UN j=0; j<outBlockMV->getNumVectors(); j++) {
400 Teuchos::ArrayRCP<SC> values = outBlockMV->getBlock(block)->getDataNonConst(j);
401 Teuchos::ArrayRCP<const SC> substractValues = substractBlockMV->getBlock(block)->getData(j);
402 if ( !vecBCType_.at(loc).compare("Dirichlet") ) {
403 for (int dd=0; dd < dofs; dd++)
404 values[ dofs * i + dd ] = substractValues[ dofs * i + dd ] - result.at(dd);
405 }
406 else if ( !vecBCType_.at(loc).compare("Dirichlet_X") ) {
407 values[ dofs * i + 0 ] = substractValues[ dofs * i + 0 ] - result.at(0);
408 }
409 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y") ) {
410 values[ dofs * i + 1 ] = substractValues[ dofs * i + 1 ] - result.at(1);
411 }
412 else if ( !vecBCType_.at(loc).compare("Dirichlet_Z") ) {
413 values[ dofs * i + 2 ] = substractValues[ dofs * i + 2 ] - result.at(2);
414 }
415 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Y") ) {
416 values[ dofs * i + 0 ] = substractValues[ dofs * i + 0 ] - result.at(0);
417 values[ dofs * i + 1 ] = substractValues[ dofs * i + 1 ] - result.at(1);
418 }
419 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Z") ) {
420 values[ dofs * i + 0 ] = substractValues[ dofs * i + 0 ] - result.at(0);
421 values[ dofs * i + 2 ] = substractValues[ dofs * i + 2 ] - result.at(2);
422 }
423 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y_Z") ) {
424 values[ dofs * i + 1 ] = substractValues[ dofs * i + 1 ] - result.at(1);
425 values[ dofs * i + 2 ] = substractValues[ dofs * i + 2 ] - result.at(2);
426 }
427 }
428 }
429 }
430 }
431 }
432 }
433 //Code below is experimental
434 for (int i=0; i<vecBCType_.size(); i++) {
435 if( vecBCType_.at(i) == "Neumann" && vecBlockID_.at(i)==block){
436 DomainPtr_Type domain = vecDomain_.at(i);
437 FEFacPtr_Type feFactory = Teuchos::rcp( new FEFac_Type() );
438 feFactory->addFE(domain);
439
440 std::vector<SC> funcParameter(3 , 0.);//0: order, 1:time, 2:surface flag (no specific flag is used in function above)
441 funcParameter[1] = t;
442 funcParameter[2] = vecFlag_.at(i);
443 int dim = domain->getDimension();
444 int dofs = vecDofs_.at(i);
445 MultiVectorPtr_Type a;
446 MultiVectorPtr_Type aUnique;
447 if (dofs>1){
448 a = Teuchos::rcp(new MultiVector_Type ( domain->getMapVecFieldRepeated() ) );
449 aUnique = Teuchos::rcp(new MultiVector_Type ( domain->getMapVecFieldUnique() ) );
450 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a, "Vector", vecBC_func_.at(i), funcParameter );
451 }
452 else{
453 a = Teuchos::rcp(new MultiVector_Type ( domain->getMapRepeated() ) );
454 aUnique = Teuchos::rcp(new MultiVector_Type ( domain->getMapUnique() ) );
455 feFactory->assemblySurfaceIntegralFlag( dim, domain->getFEType(), a, "Scalar", vecBC_func_.at(i), funcParameter );
456 }
457 aUnique->exportFromVector( a, false, "Add" );
458 MultiVectorPtr_Type mv = outBlockMV->getBlockNonConst(block);
459 mv->update(1.,*aUnique,1.);
460 MultiVectorPtr_Type mvMinus = substractBlockMV->getBlockNonConst(block);
461 mv->update(1.,*mvMinus,-1.);
462 }
463 }
464 }
465}
466
467template<class SC,class LO,class GO,class NO>
468void BCBuilder<SC,LO,GO,NO>::setAllDirichletZero(const BlockMultiVectorPtr_Type &blockMV) const{
469 vec_dbl_Type result;
470 int loc = 0;
471 for (int block = 0; block < blockMV->size(); block++) { // blocks of RHS vector
472 if(blockHasDirichletBC(block, loc)){
473 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
474 vec2D_dbl_ptr_Type vecPoints = vecDomain_.at(loc)->getPointsUnique();
475 int dim = vecDomain_.at(loc)->getDimension();
476 int dofs = vecDofs_.at(loc);
477
478 result.resize(dofs,0.);
479
480 for (int i = 0 ; i < bcFlags->size(); i++) {
481 int flag = bcFlags->at(i);
482 if(findFlag(flag,block, loc)){ // loc changed to postion corresponding to bcFlag
483 if ( !vecBCType_.at(loc).compare("Dirichlet") ||
484 !vecBCType_.at(loc).compare("Dirichlet_X") ||
485 !vecBCType_.at(loc).compare("Dirichlet_Y") ||
486 !vecBCType_.at(loc).compare("Dirichlet_Z") ||
487 !vecBCType_.at(loc).compare("Dirichlet_X_Y") ||
488 !vecBCType_.at(loc).compare("Dirichlet_X_Z") ||
489 !vecBCType_.at(loc).compare("Dirichlet_Y_Z")
490 ){
491
492 for (UN j=0; j<blockMV->getNumVectors(); j++) {
493 Teuchos::ArrayRCP<SC> values = blockMV->getBlock(block)->getDataNonConst(j);
494 if ( !vecBCType_.at(loc).compare("Dirichlet") ) {
495 for (int dd=0; dd < dofs; dd++)
496 values[ dofs * i + dd ] = result[dd];
497 }
498 else if ( !vecBCType_.at(loc).compare("Dirichlet_X") ) {
499 values[ dofs * i + 0 ] = result[0];
500 }
501 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y") ) {
502 values[ dofs * i + 1 ] = result[1];
503 }
504 else if ( !vecBCType_.at(loc).compare("Dirichlet_Z") ) {
505 values[ dofs * i + 2 ] = result[2];
506 }
507 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Y") ) {
508 values[ dofs * i + 0 ] = result[0];
509 values[ dofs * i + 1 ] = result[1];
510 }
511 else if ( !vecBCType_.at(loc).compare("Dirichlet_X_Z") ) {
512 values[ dofs * i + 0 ] = result[0];
513 values[ dofs * i + 2 ] = result[2];
514 }
515 else if ( !vecBCType_.at(loc).compare("Dirichlet_Y_Z") ) {
516 values[ dofs * i + 1 ] = result[1];
517 values[ dofs * i + 2 ] = result[2];
518 }
519 }
520 }
521 }
522 }
523 }
524 }
525}
526
527
528
529template<class SC,class LO,class GO,class NO>
530bool BCBuilder<SC,LO,GO,NO>::blockHasDirichletBC(int block) const{
531 int dummyloc;
532
533 return blockHasDirichletBC(block, dummyloc);
534}
535
536template<class SC,class LO,class GO,class NO>
537bool BCBuilder<SC,LO,GO,NO>::blockHasDirichletBC(int block, int &loc) const{
538
539 bool hasBC = false;
540 bool found = false;
541 loc = -1;
542 std::vector<int>::const_iterator it = vecBlockID_.begin();
543 while (it!=vecBlockID_.end() && !found) {
544 it = std::find(it,vecBlockID_.end(),block);
545 if (it!=vecBlockID_.end()) {
546 loc = distance(vecBlockID_.begin(),it);
547 it++;
548 if (!vecBCType_.at(loc).compare("Dirichlet") ||
549 !vecBCType_.at(loc).compare("Dirichlet_X") ||
550 !vecBCType_.at(loc).compare("Dirichlet_Y") ||
551 !vecBCType_.at(loc).compare("Dirichlet_Z") ||
552 !vecBCType_.at(loc).compare("Dirichlet_X_Y") ||
553 !vecBCType_.at(loc).compare("Dirichlet_X_Z") ||
554 !vecBCType_.at(loc).compare("Dirichlet_Y_Z")
555 ) {
556 found = true;
557 hasBC = true;
558 }
559 }
560 }
561
562 return hasBC;
563}
564
565template<class SC,class LO,class GO,class NO>
566bool BCBuilder<SC,LO,GO,NO>::findFlag(LO flag, int block, int &loc) const{
567
568#ifdef BCBuilder_TIMER
569 Teuchos::TimeMonitor FindFlagMonitor(*FindFlagTimer_);
570#endif
571
572 bool hasFlag = false;
573 bool found = false;
574 loc = -1;
575 std::vector<int>::const_iterator it = vecFlag_.begin();
576 while (it!=vecFlag_.end() && !found) {
577 it = std::find(it,vecFlag_.end(),flag);
578 if (it!=vecFlag_.end()) {
579 loc = distance(vecFlag_.begin(),it);
580 if (vecBlockID_.at(loc)==block) {
581 hasFlag = true;
582 found = true;
583 }
584 else{
585 it++;
586 }
587 }
588 }
589
590 return hasFlag;
591}
592
593template<class SC,class LO,class GO,class NO>
594void BCBuilder<SC,LO,GO,NO>::setSystem(const BlockMatrixPtr_Type &blockMatrix) const{
595
596 UN numBlocks = blockMatrix->size();
597 int loc;
598 for (UN blockRow = 0; blockRow < numBlocks; blockRow++) {
599
600#ifdef BCBuilder_TIMER
601 Teuchos::TimeMonitor SetSystemRowMonitor(*SetSystemRowTimer_);
602#endif
603 bool boolBlockHasDirichlet = false;
604
605 {
606#ifdef BCBuilder_TIMER
607 Teuchos::TimeMonitor BlockRowHasDirichletMonitor(*BlockRowHasDirichletTimer_);
608#endif
609
610 boolBlockHasDirichlet = blockHasDirichletBC(blockRow,loc);
611 }
612
613 if (boolBlockHasDirichlet){
614 for (UN blockCol = 0; blockCol < numBlocks ; blockCol++) {
615 if (blockMatrix->blockExists(blockRow, blockCol)) {
616 MatrixPtr_Type matrix = blockMatrix->getBlock(blockRow, blockCol);
617 setDirichletBC(matrix, loc, blockRow, blockRow == blockCol);
618 } else if (blockRow == blockCol) {
619 // If the current block is a diagonal block but does not exist we need to create it and insert the
620 // value 1 in all diagonal entries of the block corresponding to a Dirichlet boundary condition. We
621 // only need max. one entry per row.
622
623 // To get the correct domain from vecDomain_ we need to know the index of any (in this case the
624 // first) entry that was done with addBC for this block
625 const auto it = std::find(this->vecBlockID_.begin(), this->vecBlockID_.end(), blockRow);
626 std::cout << " loc " << loc << " it - something " << it - this->vecBlockID_.begin() << std::endl;
627 auto matrixMap = this->vecDomain_.at(loc)->getMapUnique();
628 // // Use Xpetra::MatrixFactory to build a matrix with known row and column maps
629 auto tpetraMatrix = Teuchos::rcp(new Tpetra::CrsMatrix<SC,LO,GO,NO>(matrixMap->getTpetraMap(), matrixMap->getTpetraMap(), 1)); // Tpetra::MatrixFactory<SC, LO, GO, NO>::Build(matrixMap->getTpetraMap(), matrixMap->getTpetraMap(), 1);
630
631 MatrixPtr_Type matrix = Teuchos::rcp(new Matrix_Type(tpetraMatrix));
632
633 // Fill the matrix with zeros on the diagonal
634 Teuchos::Array<LO> colIndex(1, 0);
635 Teuchos::Array<SC> val(1, Teuchos::ScalarTraits<SC>::zero());
636 for (auto i = 0; i < matrixMap->getNodeNumElements(); i++) {
637 colIndex[0] = i;
638 // Since the matrix has a column map we can use insertLocalValues(). This is more efficient!
639 matrix->insertLocalValues(i, colIndex, val);
640 }
641 matrix->fillComplete(matrixMap, matrixMap);
642 setDirichletBC(matrix, loc, blockRow, true);
643 blockMatrix->addBlock(matrix, blockRow, blockCol);
644 }
645 }
646 }
647 }
648}
649
650template<class SC,class LO,class GO,class NO>
651void BCBuilder<SC,LO,GO,NO>::setDirichletBC(const MatrixPtr_Type &matrix, int loc, int blockRow, bool isDiagonalBlock) const{
652
653 matrix->resumeFill();
654 bool isDirichlet;
655 UN dofsPerNode = vecDofs_.at(loc);
656 vec_int_ptr_Type bcFlags = vecDomain_.at(loc)->getBCFlagUnique();
657
658 for (LO i=0; i<bcFlags->size(); i++) { // flags are for nodes.
659 int flag = bcFlags->at(i);
660 isDirichlet = false;
661 int locThisFlag;
662
663 if( findFlag(flag, blockRow, locThisFlag) ){
664
665 if( !vecBCType_.at(locThisFlag).compare("Dirichlet") ||
666 !vecBCType_.at(locThisFlag).compare("Dirichlet_X") ||
667 !vecBCType_.at(locThisFlag).compare("Dirichlet_Y") ||
668 !vecBCType_.at(locThisFlag).compare("Dirichlet_Z") ||
669 !vecBCType_.at(locThisFlag).compare("Dirichlet_X_Y") ||
670 !vecBCType_.at(locThisFlag).compare("Dirichlet_X_Z") ||
671 !vecBCType_.at(locThisFlag).compare("Dirichlet_Y_Z") )
672 isDirichlet = true;
673 if (isDirichlet) {
674
675 if (isDiagonalBlock)
676 setLocalRowOne(matrix, i, dofsPerNode, locThisFlag );
677 else
678 setLocalRowZero(matrix, i, dofsPerNode, locThisFlag );
679 }
680 }
681 }
682 matrix->fillComplete( matrix->getMap("domain"), matrix->getMap("range") );
683}
684
685template<class SC,class LO,class GO,class NO>
686void BCBuilder<SC,LO,GO,NO>::setLocalRowOne(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc) const{
687
688 Teuchos::ArrayView<const SC> valuesOld;
689 Teuchos::ArrayView<const LO> indices;
690
691 MapConstPtr_Type colMap = matrix->getMap("col");
692 LO localDof = (LO) (dofsPerNode * localNode);
693 for (UN dof=0; dof<dofsPerNode; dof++) {
694 if ( vecBCType_.at(loc) == "Dirichlet" ||
695 (vecBCType_.at(loc) == "Dirichlet_X" && dof==0 ) ||
696 (vecBCType_.at(loc) == "Dirichlet_Y" && dof==1 ) ||
697 (vecBCType_.at(loc) == "Dirichlet_Z" && dof==2 ) ||
698 (vecBCType_.at(loc) == "Dirichlet_X_Y" && dof!=2 )||
699 (vecBCType_.at(loc) == "Dirichlet_X_Z" && dof!=1 )||
700 (vecBCType_.at(loc) == "Dirichlet_Y_Z" && dof!=0 )
701 ) {
702 // std::cout << " Setting Dirichlet Row 1 for node " << localDof << " of type " << vecBCType_.at(loc) <<endl;
703 GO globalDof = matrix->getMap()->getGlobalElement( localDof );
704 matrix->getLocalRowView(localDof, indices, valuesOld);
705 Teuchos::Array<SC> values( valuesOld.size(), Teuchos::ScalarTraits<SC>::zero() );
706 bool setOne = false;
707 for (UN j=0; j<indices.size() && !setOne; j++) {
708 if ( colMap->getGlobalElement( indices[j] ) == globalDof ){
709 values[j] = Teuchos::ScalarTraits<SC>::one();
710 setOne = true;
711 }
712 }
713 matrix->replaceLocalValues(localDof, indices(), values());
714 }
715 localDof++;
716 }
717
718}
719
720template<class SC,class LO,class GO,class NO>
721void BCBuilder<SC,LO,GO,NO>::setLocalRowZero(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc) const{
722
723 Teuchos::ArrayView<const SC> valuesOld;
724 Teuchos::ArrayView<const LO> indices;
725 LO localDof = (LO) (dofsPerNode * localNode);
726 for (UN dof=0; dof<dofsPerNode; dof++) {
727 if ( vecBCType_.at(loc) == "Dirichlet" ||
728 (vecBCType_.at(loc) == "Dirichlet_X" && dof==0 ) ||
729 (vecBCType_.at(loc) == "Dirichlet_Y" && dof==1 ) ||
730 (vecBCType_.at(loc) == "Dirichlet_Z" && dof==2 ) ||
731 (vecBCType_.at(loc) == "Dirichlet_X_Y" && dof!=2 )||
732 (vecBCType_.at(loc) == "Dirichlet_X_Z" && dof!=1 )||
733 (vecBCType_.at(loc) == "Dirichlet_Y_Z" && dof!=0 )
734 ) {
735 matrix->getLocalRowView(localDof, indices, valuesOld);
736 Teuchos::Array<SC> values( valuesOld.size(), Teuchos::ScalarTraits<SC>::zero() );
737 matrix->replaceLocalValues(localDof, indices(), values());
738 }
739 localDof++;
740 }
741}
742
743
744
745template<class SC,class LO,class GO,class NO>
746int BCBuilder<SC,LO,GO,NO>::dofsPerNodeAtBlock(int block) {
747 std::vector<int>::iterator it;
748 it = find(vecBlockID_.begin(),vecBlockID_.end(),block);
749#ifdef ASSERTS_WARNINGS
750 MYASSERT(it!=vecBlockID_.end(),"Requested dof information is not known to BCBuilder class.");
751#endif
752
753 return vecDofs_.at(distance(vecBlockID_.begin(),it));
754}
755}
756
757#endif
FEDD::BCBuilder::setBCMinusVector
void setBCMinusVector(const BlockMultiVectorPtr_Type &outBlockMV, const BlockMultiVectorPtr_Type &substractBlockMV, double t=0.) const
Definition BCBuilder_def.hpp:238
FEDD::BCBuilder::setVectorMinusBC
void setVectorMinusBC(const BlockMultiVectorPtr_Type &outBlockMV, const BlockMultiVectorPtr_Type &substractBlockMV, double t=0.) const
Definition BCBuilder_def.hpp:354
FEDD::BCBuilder::setRHS
void setRHS(const BlockMultiVectorPtr_Type &blockMV, double t=0.) const
Setting boundary conditions to (block)vector.
Definition BCBuilder_def.hpp:97
FEDD::BCBuilder::set
void set(const BlockMatrixPtr_Type &blockMatrix, const BlockMultiVectorPtr_Type &b, double t=0.) const
Setting bundary condtions to problem.
Definition BCBuilder_def.hpp:88
FEDD::BCBuilder::setSystem
void setSystem(const BlockMatrixPtr_Type &blockMatrix) const
Set boundary conditions to system.
Definition BCBuilder_def.hpp:594
FEDD::BCBuilder::addBC
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:42
FEDD::BCBuilder::setLocalRowZero
void setLocalRowZero(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc) const
Definition BCBuilder_def.hpp:721
FEDD::BCBuilder< SC, LO, GO, NO >::setDirichletBoundaryFromExternal
void setDirichletBoundaryFromExternal(Teuchos::ArrayRCP< SC > &values, LO index, int loc, double time, std::string type, Teuchos::ArrayRCP< SC > valuesSubstract=Teuchos::null) const
FEDD::BCBuilder::setDirichletBC
void setDirichletBC(const MatrixPtr_Type &matrix, int loc, int blockRow, bool isDiagonalBlock) const
Definition BCBuilder_def.hpp:651
FEDD::BCBuilder::setAllDirichletZero
void setAllDirichletZero(const BlockMultiVectorPtr_Type &blockMV) const
Definition BCBuilder_def.hpp:468
FEDD::BCBuilder< SC, LO, GO, NO >::blockHasDirichletBC
bool blockHasDirichletBC(int block) const
FEDD::BCBuilder< SC, LO, GO, NO >::findFlag
bool findFlag(LO flag, int block, int &loc) const
FEDD::BCBuilder::setLocalRowOne
void setLocalRowOne(const MatrixPtr_Type &matrix, LO localNode, UN dofsPerNode, int loc) const
Definition BCBuilder_def.hpp:686
FEDD::BCBuilder::dofsPerNodeAtBlock
int dofsPerNodeAtBlock(int block)
Definition BCBuilder_def.hpp:746
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5