00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <BoundaryConditionDiscretizationSpectralConform.hpp>
00021 #include <ErrorHandler.hpp>
00022
00023 #include <VariationalProblem.hpp>
00024 #include <DegreeOfFreedomSet.hpp>
00025 #include <VariationalBorderOperatorAlphaUV.hpp>
00026 #include <VariationalBorderOperatorFV.hpp>
00027 #include <VariationalLinearOperator.hpp>
00028
00029 #include <BaseMatrix.hpp>
00030 #include <BaseVector.hpp>
00031
00032 #include <Discretization.hpp>
00033 #include <ScalarDiscretizationTypeSpectral.hpp>
00034
00035 #include <BoundaryReferences.hpp>
00036
00037 #include <UnAssembledMatrix.hpp>
00038
00039 #include <DoubleHashedMatrix.hpp>
00040
00041
00042 void BoundaryConditionDiscretizationSpectralConform::
00043 getDiagonal(BaseVector& X) const
00044 {
00045 Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(X);
00046
00047 switch(this->problem().type()) {
00048 case Problem::variational: {
00049 const VariationalProblem& variationalProblem
00050 = dynamic_cast<const VariationalProblem&>(this->problem());
00051
00052 for (VariationalProblem::bilinearBorderOperatorConst_iterator
00053 iBorderOperator = variationalProblem.beginBilinearBorderOperator();
00054 iBorderOperator != variationalProblem.endBilinearBorderOperator();
00055 ++iBorderOperator) {
00056
00057 const VariationalBilinearBorderOperator& borderOperator = **iBorderOperator;
00058 switch (borderOperator.type()) {
00059 case VariationalBilinearBorderOperator::alphaUV: {
00060 const VariationalBorderOperatorAlphaUV& operatorAlphaUV
00061 = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(borderOperator);
00062 ConstReferenceCounting<Boundary> boundary = operatorAlphaUV.boundary();
00063
00064 switch(boundary->type()) {
00065 case Boundary::references: {
00066 const BoundaryReferences& boundaryReference =
00067 dynamic_cast<const BoundaryReferences&>(*boundary);
00068
00069 for (BoundaryReferences::ReferencesSet::const_iterator
00070 iBoundary = boundaryReference.references().begin();
00071 iBoundary != boundaryReference.references().end(); ++iBoundary) {
00072
00073 const ScalarFunctionBase& alpha = operatorAlphaUV.alpha();
00074 const size_t boundaryNumber = *iBoundary;
00075
00076 const size_t& numberOfUnknown = this->problem().numberOfUnknown();
00077
00078 Vector<TinyVector<3,Vector<Vector<real_t> > > > quadratureBaseValuesU(numberOfUnknown);
00079 TinyVector<3,size_t> direction;
00080 for (size_t i=0; i<3; ++i) {
00081 direction[i] = (boundaryNumber/2+i+1) % 3;
00082 }
00083
00084 for (size_t unknownNumber =0; unknownNumber < numberOfUnknown;++ unknownNumber){
00085
00086 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
00087 = __basisU[unknownNumber];
00088
00089 for (size_t i=0; i<2; ++i) {
00090 quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00091 for (size_t j=0; j < (*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00092
00093 real_t xj=(*__transformU[unknownNumber][direction[i]]).inverse((*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j)));
00094
00095 quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[direction[i]]->dimension());
00096
00097 unknownBasis[direction[i]]->getValues(xj,quadratureBaseValuesU[unknownNumber][i][j]);
00098 }
00099 }
00100
00101 quadratureBaseValuesU[unknownNumber][2].resize(2);
00102 for (size_t j=0; j<2 ; ++j) {
00103 real_t xj= -1. + 2.*j;
00104 quadratureBaseValuesU[unknownNumber][2][j].resize(unknownBasis[direction[2]]->dimension());
00105 unknownBasis[direction[2]]->getValues(xj,quadratureBaseValuesU[unknownNumber][2][j]);
00106 }
00107 }
00108
00109 TinyVector<2, Vector<real_t> >nodes;
00110 for (size_t i=0; i<2; ++i) {
00111 nodes[i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00112 for (size_t j=0; j<(*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00113 nodes[i][j] =(*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j));
00114 }
00115 }
00116 const real_t jacobianDet = __transform[direction[0]]->determinant()*
00117 __transform[direction[1]]->determinant();
00118
00119 Vector<Vector<real_t> > alpha_rs(__gaussLobatto[direction[0]]->numberOfPoints());
00120 for(size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00121 alpha_rs[r].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00122 for(size_t s=0; s< __gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00123 alpha_rs[r][s]=0;
00124 }
00125 }
00126
00127 for (size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00128 const real_t& x = nodes[0][r];
00129 for (size_t s=0; s <__gaussLobatto[direction[1]]->numberOfPoints(); ++s) {
00130 const real_t& y = nodes[1][s];
00131
00132 switch(boundaryNumber) {
00133 case 0: {
00134 alpha_rs[r][s] = alpha(__mesh.shape().a(0), x, y);
00135 break;
00136 }
00137 case 1: {
00138 alpha_rs[r][s] = alpha(__mesh.shape().b(0), x, y);
00139 break;
00140 }
00141 case 2: {
00142 alpha_rs[r][s] = alpha(y, __mesh.shape().a(1), x);
00143 break;
00144 }
00145 case 3: {
00146 alpha_rs[r][s] = alpha(y, __mesh.shape().b(1), x);
00147 break;
00148 }
00149 case 4: {
00150 alpha_rs[r][s] = alpha(x, y, __mesh.shape().a(2));
00151 break;
00152 }
00153 case 5: {
00154 alpha_rs[r][s] = alpha(x, y, __mesh.shape().b(2));
00155 break;
00156 }
00157 }
00158 }
00159 }
00160
00161
00162 const TinyVector<3,size_t>& dofShape = __dofShape[operatorAlphaUV.unknownNumber()];
00163
00164 TinyVector<3,size_t> I;
00165
00166 for(size_t unknownNumber=0;unknownNumber< numberOfUnknown; unknownNumber++) {
00167
00168 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
00169 = __basisU[unknownNumber];
00170
00171 for(size_t m=0; m< unknownBasis[direction[2]]->dimension(); ++m){
00172 I[direction[2]]=m;
00173 for(size_t n=0; n< unknownBasis[direction[0]]->dimension(); ++n){
00174 I[direction[0]]=n;
00175 for(size_t p=0; p<unknownBasis[direction[1]]->dimension(); ++p){
00176 I[direction[1]]=p;
00177 size_t dofNumber = I[0]*dofShape[1]*dofShape[2] + I[1]*dofShape[2] + I[2];
00178
00179 for(size_t r=0; r< __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00180 for(size_t s=0; s < __gaussLobatto[direction[1]]->numberOfPoints(); ++s) {
00181 v[__degreeOfFreedomSet(operatorAlphaUV.testFunctionNumber(),dofNumber)]
00182 += quadratureBaseValuesU[operatorAlphaUV.unknownNumber()][2][boundaryNumber%2][m]
00183 *quadratureBaseValuesU[operatorAlphaUV.testFunctionNumber()][2][boundaryNumber%2][m]
00184 *quadratureBaseValuesU[operatorAlphaUV.unknownNumber()][0][n]
00185 *quadratureBaseValuesU[operatorAlphaUV.testFunctionNumber()][0][n]
00186 *quadratureBaseValuesU[operatorAlphaUV.unknownNumber()][1][p]
00187 *quadratureBaseValuesU[operatorAlphaUV.testFunctionNumber()][1][p]
00188 *__gaussLobatto[direction[0]]->weight(r)
00189 *__gaussLobatto[direction[1]]->weight(s)
00190 *alpha_rs[r][s]*jacobianDet;
00191 }
00192 }
00193 }
00194 }
00195 }
00196 }
00197 }
00198
00199 break;
00200 }
00201 default: {
00202 throw ErrorHandler(__FILE__,__LINE__,
00203 "not implemented",
00204 ErrorHandler::unexpected);
00205 }
00206 }
00207 break;
00208 }
00209 default: {
00210 throw ErrorHandler(__FILE__,__LINE__,
00211 "not implemented",
00212 ErrorHandler::unexpected);
00213 }
00214 }
00215 }
00216 break;
00217 }
00218 case Problem::pde: {
00219 break;
00220 }
00221 default: {
00222 throw ErrorHandler(__FILE__,__LINE__,
00223 "not implemented",
00224 ErrorHandler::unexpected);
00225 }
00226 }
00227 }
00228
00229
00230 void BoundaryConditionDiscretizationSpectralConform::
00231 setMatrix(ReferenceCounting<BaseMatrix> givenA,
00232 ReferenceCounting<BaseVector> b) const
00233 {
00234 switch((*givenA).type()) {
00235 case BaseMatrix::doubleHashedMatrix: {
00236 throw ErrorHandler(__FILE__,__LINE__,
00237 "not implemented",
00238 ErrorHandler::unexpected);
00239 break;
00240 }
00241 case BaseMatrix::unAssembled: {
00242 UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(*givenA);
00243 A.setBoundaryConditions(this);
00244 break;
00245 }
00246
00247 default: {
00248 throw ErrorHandler(__FILE__,__LINE__,
00249 "unexected matrix type",
00250 ErrorHandler::unexpected);
00251 }
00252 }
00253 }
00254
00255
00256 void BoundaryConditionDiscretizationSpectralConform::
00257 setSecondMember(ReferenceCounting<BaseMatrix> givenA,
00258 ReferenceCounting<BaseVector> givenb) const
00259 {
00260 ffout(2)<<"- assembling second membre-boundary\n";
00261 Vector<real_t>& b= (static_cast<Vector<real_t>&>(*givenb));
00262
00263 switch(this->problem().type()) {
00264 case Problem::variational: {
00265 const VariationalProblem& variationalProblem
00266 = dynamic_cast<const VariationalProblem&>(this->problem());
00267
00268 for (VariationalProblem::linearBorderOperatorConst_iterator
00269 iLinearBorderOperator = variationalProblem.beginLinearBorderOperator();
00270 iLinearBorderOperator != variationalProblem.endLinearBorderOperator();
00271 ++iLinearBorderOperator) {
00272
00273 const VariationalLinearBorderOperator& linearBorderOperator = **iLinearBorderOperator;
00274 switch (linearBorderOperator.type()) {
00275 case VariationalLinearBorderOperator::FV: {
00276 const VariationalBorderOperatorFV& operatorFV
00277 = dynamic_cast<const VariationalBorderOperatorFV&>(linearBorderOperator);
00278 ConstReferenceCounting<Boundary> boundary = operatorFV.boundary();
00279
00280 switch(boundary->type()) {
00281 case Boundary::references: {
00282 const BoundaryReferences& boundaryReference =
00283 dynamic_cast<const BoundaryReferences&>(*boundary);
00284
00285 for (BoundaryReferences::ReferencesSet::const_iterator
00286 iBoundary = boundaryReference.references().begin();
00287 iBoundary != boundaryReference.references().end(); ++iBoundary) {
00288
00289 const ScalarFunctionBase& g = operatorFV.f();
00290 const size_t boundaryNumber = *iBoundary;
00291
00292 const size_t& numberOfUnknown = this->problem().numberOfUnknown();
00293
00294 TinyVector<3, size_t> direction;
00295 for (size_t i=0; i<3; ++i) {
00296 direction[i] = (boundaryNumber/2+i+1) % 3;
00297 }
00298
00299 Vector<TinyVector<3,Vector<Vector<real_t> > > > quadratureBaseValues(numberOfUnknown);
00300 for (size_t testFunctionNumber =0; testFunctionNumber < numberOfUnknown; ++ testFunctionNumber){
00301 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& testBasis
00302 = __basisU[testFunctionNumber];
00303
00304 for (size_t i=0; i<2; ++i) {
00305 quadratureBaseValues[testFunctionNumber][i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00306
00307 for (size_t j=0; j < (*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00308 real_t xj= (*__transformU[testFunctionNumber][direction[i]]).inverse((*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j)));
00309
00310 quadratureBaseValues[testFunctionNumber][i][j].resize(testBasis[direction[i]]->dimension());
00311 testBasis[direction[i]]->getValues(xj,quadratureBaseValues[testFunctionNumber][i][j]);
00312 }
00313 }
00314
00315 quadratureBaseValues[testFunctionNumber][2].resize(2);
00316 for (size_t j=0; j<2 ; ++j) {
00317 real_t xj= -1.+2.*j;
00318 quadratureBaseValues[testFunctionNumber][2][j].resize(testBasis[direction[2]]->dimension());
00319 testBasis[direction[2]]->getValues(xj,quadratureBaseValues[testFunctionNumber][2][j]);
00320 }
00321 }
00322
00323
00324 TinyVector<2, Vector<real_t> >nodes;
00325 for (size_t i=0; i<2; ++i) {
00326 nodes[i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00327 for (size_t j=0; j< (*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00328 nodes[i][j] =(*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j));
00329 }
00330 }
00331
00332 const real_t jacobianDet = (*__transform[direction[0]]).determinant()*
00333 (*__transform[direction[1]]).determinant();
00334
00335 for (size_t i=0; i < (*__gaussLobatto[direction[0]]).numberOfPoints(); ++i) {
00336 const real_t& x = nodes[0][i];
00337 const real_t& wi = (*__gaussLobatto[direction[0]]).weight(i)
00338 ;
00339
00340 for (size_t j=0; j < (*__gaussLobatto[direction[1]]).numberOfPoints(); ++j) {
00341 const real_t& y = nodes[1][j];
00342 const real_t& wij = (*__gaussLobatto[direction[1]]).weight(j)*wi
00343 ;
00344
00345 real_t gValue = wij* jacobianDet;
00346 switch(boundaryNumber) {
00347 case 0: {
00348 gValue *= g(__mesh.shape().a(0), x, y);
00349 break;
00350 }
00351 case 1: {
00352 gValue *= g(__mesh.shape().b(0), x, y);
00353 break;
00354 }
00355 case 2: {
00356 gValue *= g(y, __mesh.shape().a(1), x);
00357 break;
00358 }
00359 case 3: {
00360 gValue *= g(y, __mesh.shape().b(1), x);
00361 break;
00362 }
00363 case 4: {
00364 gValue *= g(x, y, __mesh.shape().a(2));
00365 break;
00366 }
00367 case 5: {
00368 gValue *= g(x, y, __mesh.shape().b(2));
00369 break;
00370 }
00371 default: {
00372 throw ErrorHandler(__FILE__,__LINE__,
00373 "not implemented",
00374 ErrorHandler::unexpected);
00375 }
00376 }
00377
00378 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& testBasis
00379 = __basisU[operatorFV.testFunctionNumber()];
00380
00381 const TinyVector<3,size_t>& dofShape =
00382 __dofShape[operatorFV.testFunctionNumber()];
00383
00384 TinyVector<3,size_t> J;
00385 for (size_t ml=0; ml < testBasis[direction[0]]->dimension(); ++ml){
00386 J[direction[0]] = ml;
00387 for (size_t pl=0; pl < testBasis[direction[1]]->dimension(); ++pl){
00388 J[direction[1]] = pl;
00389 for (size_t ql=0; ql < testBasis[direction[2]]->dimension(); ++ql){
00390 J[direction[2]] = ql;
00391 const size_t l =
00392 J[0]*dofShape[1]*dofShape[2] + J[1]*dofShape[2] + J[2];
00393 b[__degreeOfFreedomSet(operatorFV.testFunctionNumber(),l)]
00394
00395 +=quadratureBaseValues[operatorFV.testFunctionNumber()][0][i][ml]
00396 * quadratureBaseValues[operatorFV.testFunctionNumber()][1][j][pl]
00397 * quadratureBaseValues[operatorFV.testFunctionNumber()][2][boundaryNumber%2][ql]
00398 *gValue
00399 ;
00400 }
00401 }
00402 }
00403 }
00404 }
00405 }
00406
00407 break;
00408 }
00409 default: {
00410 throw ErrorHandler(__FILE__,__LINE__,
00411 "not implemented",
00412 ErrorHandler::unexpected);
00413 }
00414
00415 }
00416 break;
00417 }
00418 default: {
00419 throw ErrorHandler(__FILE__,__LINE__,
00420 "not implemented",
00421 ErrorHandler::unexpected);
00422 }
00423 }
00424 }
00425
00426 break;
00427 }
00428
00429 default: {
00430 throw ErrorHandler(__FILE__,__LINE__,
00431 "not implemented",
00432 ErrorHandler::unexpected);
00433 }
00434 }
00435 }
00436
00437
00438 void BoundaryConditionDiscretizationSpectralConform::
00439 timesX(const BaseVector& X, BaseVector& Z) const
00440 {
00441 const Vector<real_t>& u = dynamic_cast<const Vector<real_t>&>(X);
00442 Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(Z);
00443
00444 switch(this->problem().type()) {
00445 case Problem::variational: {
00446 const VariationalProblem& variationalProblem
00447 = dynamic_cast<const VariationalProblem&>(this->problem());
00448
00449 for (VariationalProblem::bilinearBorderOperatorConst_iterator
00450 iBorderOperator = variationalProblem.beginBilinearBorderOperator();
00451 iBorderOperator != variationalProblem.endBilinearBorderOperator();
00452 ++iBorderOperator) {
00453
00454 const VariationalBilinearBorderOperator& borderOperator = **iBorderOperator;
00455 switch (borderOperator.type()) {
00456 case VariationalBilinearBorderOperator::alphaUV: {
00457
00458 const VariationalBorderOperatorAlphaUV& operatorAlphaUV
00459 = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(borderOperator);
00460 ConstReferenceCounting<Boundary> boundary = operatorAlphaUV.boundary();
00461
00462 switch(boundary->type()) {
00463 case Boundary::references: {
00464 const BoundaryReferences& boundaryReference = dynamic_cast<const BoundaryReferences&>(*boundary);
00465 for (BoundaryReferences::ReferencesSet::const_iterator
00466 iBoundary = boundaryReference.references().begin();
00467 iBoundary != boundaryReference.references().end(); ++iBoundary) {
00468
00469 const ScalarFunctionBase& alpha = operatorAlphaUV.alpha();
00470 const size_t boundaryNumber = *iBoundary;
00471
00472 const size_t& numberOfUnknown = this->problem().numberOfUnknown();
00473
00474 Vector<TinyVector<3,Vector<Vector<real_t> > > > quadratureBaseValuesU(numberOfUnknown);
00475 TinyVector<3,size_t> direction;
00476 for (size_t i=0; i<3; ++i) {
00477 direction[i] = (boundaryNumber/2+i+1) % 3;
00478 }
00479
00480 for (size_t unknownNumber =0; unknownNumber < numberOfUnknown;++ unknownNumber){
00481
00482 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >& unknownBasis
00483 = __basisU[unknownNumber];
00484
00485 for (size_t i=0; i<2; ++i) {
00486 quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00487 for (size_t j=0; j < (*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00488
00489 real_t xj=(*__transformU[unknownNumber][direction[i]]).inverse((*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j)));
00490
00491 quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[direction[i]]->dimension());
00492
00493 unknownBasis[direction[i]]->getValues(xj,quadratureBaseValuesU[unknownNumber][i][j]);
00494 }
00495 }
00496
00497 quadratureBaseValuesU[unknownNumber][2].resize(2);
00498 for (size_t j=0; j<2 ; ++j) {
00499 real_t xj= -1. + 2.*j;
00500 quadratureBaseValuesU[unknownNumber][2][j].resize(unknownBasis[direction[2]]->dimension());
00501 unknownBasis[direction[2]]->getValues(xj,quadratureBaseValuesU[unknownNumber][2][j]);
00502 }
00503 }
00504
00505
00506
00507 TinyVector<2, Vector<real_t> >nodes;
00508 for (size_t i=0; i<2; ++i) {
00509 nodes[i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00510 for (size_t j=0; j<(*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00511 nodes[i][j] =(*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j));
00512 }
00513 }
00514 const real_t jacobianDet = __transform[direction[0]]->determinant()*
00515 __transform[direction[1]]->determinant();
00516
00517 const size_t unknownNumber= operatorAlphaUV.unknownNumber();
00518 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >&unknownBasis
00519 = __basisU[unknownNumber];
00520
00521 Vector<Vector<real_t> > u_jk(unknownBasis[direction[0]]->dimension());
00522 for(size_t j=0; j< unknownBasis[direction[0]]->dimension(); ++j){
00523 u_jk[j].resize(unknownBasis[direction[1]]->dimension());
00524 for(size_t k=0; k< unknownBasis[direction[1]]->dimension(); ++k){
00525 u_jk[j][k]=0;
00526 }
00527 }
00528
00529 {
00530 const TinyVector<3,size_t>& dofShapeU = __dofShape[unknownNumber];
00531 TinyVector<3,size_t> J;
00532 for(size_t j=0; j< unknownBasis[direction[0]]->dimension(); ++j){
00533 J[direction[0]]=j;
00534 for(size_t k=0; k< unknownBasis[direction[1]]->dimension(); ++k){
00535 J[direction[1]]=k;
00536 real_t& u_jk_value = u_jk[j][k];
00537 for(size_t i=0; i< unknownBasis[direction[2]]->dimension(); ++i){
00538 J[direction[2]]=i;
00539 const size_t dofNumberU =
00540 J[0]*dofShapeU[1]*dofShapeU[2] + J[1] *dofShapeU[2] +J[2];
00541 u_jk_value
00542 += quadratureBaseValuesU[unknownNumber][2][boundaryNumber%2][i]
00543 *u[__degreeOfFreedomSet(operatorAlphaUV.unknownNumber(),dofNumberU)];
00544 }
00545 }
00546 }
00547 }
00548
00549 Vector<Vector<real_t> > u_js(unknownBasis[direction[0]]->dimension());
00550 for(size_t j=0; j<unknownBasis[direction[0]]->dimension(); ++j){
00551 u_js[j].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00552 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00553 u_js[j][s]=0;
00554 }
00555 }
00556
00557 for(size_t j=0; j< unknownBasis[direction[0]]->dimension(); ++j){
00558 for(size_t s=0; s< __gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00559 real_t& u_js_value = u_js[j][s];
00560 for(size_t k=0; k< unknownBasis[direction[1]]->dimension(); ++k){
00561 u_js_value
00562 += quadratureBaseValuesU[unknownNumber][1][s][k]* u_jk[j][k];
00563 }
00564 }
00565 }
00566
00567 Vector<Vector<real_t> > u_rs(__gaussLobatto[direction[0]]->numberOfPoints());
00568 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00569 u_rs[r].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00570 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00571 u_rs[r][s]=0;
00572 }
00573 }
00574
00575 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00576 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00577 real_t& u_rs_value = u_rs[r][s];
00578 for(size_t j=0; j<unknownBasis[direction[0]]->dimension(); ++j){
00579 u_rs_value
00580 += quadratureBaseValuesU[unknownNumber][0][r][j]*u_js[j][s];
00581 }
00582 }
00583 }
00584
00585
00586 Vector<Vector<real_t> > alpha_rs(__gaussLobatto[direction[0]]->numberOfPoints());
00587 for(size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00588 alpha_rs[r].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00589 for(size_t s=0; s< __gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00590 alpha_rs[r][s]=0;
00591 }
00592 }
00593
00594 for (size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00595 const real_t& x = nodes[0][r];
00596 for (size_t s=0; s <__gaussLobatto[direction[1]]->numberOfPoints(); ++s) {
00597 const real_t& y = nodes[1][s];
00598
00599
00600 switch(boundaryNumber) {
00601 case 0: {
00602 alpha_rs[r][s] = alpha(__mesh.shape().a(0), x, y);
00603 break;
00604 }
00605 case 1: {
00606 alpha_rs[r][s] = alpha(__mesh.shape().b(0), x, y);
00607 break;
00608 }
00609 case 2: {
00610 alpha_rs[r][s] = alpha(y, __mesh.shape().a(1), x);
00611 break;
00612 }
00613 case 3: {
00614 alpha_rs[r][s] = alpha(y, __mesh.shape().b(1), x);
00615 break;
00616 }
00617 case 4: {
00618 alpha_rs[r][s] = alpha(x, y, __mesh.shape().a(2));
00619 break;
00620 }
00621 case 5: {
00622 alpha_rs[r][s] = alpha(x, y, __mesh.shape().b(2));
00623 break;
00624 }
00625 }
00626 }
00627 }
00628
00629 const size_t testFunctionNumber= operatorAlphaUV.testFunctionNumber();
00630 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >&testBasis
00631 = __basisU[testFunctionNumber];
00632 Vector<Vector<real_t> > v_rp(__gaussLobatto[direction[0]]->numberOfPoints());
00633 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00634 v_rp[r].resize(testBasis[direction[1]]->dimension());
00635 for(size_t p=0; p<testBasis[direction[1]]->dimension(); ++p){
00636 v_rp[r][p]=0;
00637 }
00638 }
00639
00640 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00641 for(size_t p=0; p<testBasis[direction[1]]->dimension(); ++p){
00642 real_t& v_rp_value = v_rp[r][p];
00643 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00644 v_rp_value
00645 += quadratureBaseValuesU[testFunctionNumber][1][s][p]
00646 *u_rs[r][s]
00647 *__gaussLobatto[direction[1]]->weight(s)
00648 * alpha_rs[r][s]*jacobianDet;
00649 }
00650 }
00651 }
00652
00653 Vector<Vector<real_t> > v_np(testBasis[direction[0]]->dimension());
00654 for(size_t n=0; n<testBasis[direction[0]]->dimension(); ++n){
00655 v_np[n].resize(testBasis[direction[1]]->dimension());
00656 for(size_t p=0; p<testBasis[direction[1]]->dimension(); ++p){
00657 v_np[n][p]=0;
00658 }
00659 }
00660
00661 for(size_t n=0; n<testBasis[direction[0]]->dimension(); ++n){
00662 for(size_t p=0; p<testBasis[direction[1]]->dimension(); ++p){
00663 real_t& v_np_value = v_np[n][p];
00664 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00665 v_np_value
00666 += quadratureBaseValuesU[testFunctionNumber][0][r][n]
00667 *v_rp[r][p]*
00668 __gaussLobatto[direction[0]]->weight(r);
00669 }
00670 }
00671 }
00672
00673 const TinyVector<3,size_t>& dofShapeV = __dofShape[testFunctionNumber];
00674 TinyVector<3,size_t> I;
00675
00676 for(size_t n=0; n<testBasis[direction[0]]->dimension(); ++n){
00677 I[direction[0]]=n;
00678 for(size_t p=0; p<testBasis[direction[1]]->dimension(); ++p){
00679 I[direction[1]]=p;
00680 for(size_t m=0; m<testBasis[direction[2]]->dimension(); ++m){
00681 I[direction[2]]=m;
00682 size_t dofNumberV =
00683 I[0]*dofShapeV[1]*dofShapeV[2] + I[1]*dofShapeV[2] +I[2];
00684
00685 v[__degreeOfFreedomSet(operatorAlphaUV.testFunctionNumber(),dofNumberV)]
00686 += quadratureBaseValuesU[testFunctionNumber][2][boundaryNumber%2][m]
00687 *v_np[n][p];
00688 }
00689 }
00690 }
00691 }
00692 break;
00693 }
00694 default: {
00695 throw ErrorHandler(__FILE__,__LINE__,
00696 "not implemented",
00697 ErrorHandler::unexpected);
00698 }
00699 }
00700 break;
00701 }
00702 default: {
00703 throw ErrorHandler(__FILE__,__LINE__,
00704 "not implemented",
00705 ErrorHandler::unexpected);
00706 }
00707 }
00708 }
00709 break;
00710 }
00711 case Problem::pde: {
00712 break;
00713 }
00714 default: {
00715 throw ErrorHandler(__FILE__,__LINE__,
00716 "not implemented",
00717 ErrorHandler::unexpected);
00718 }
00719 }
00720 }
00721
00722 void BoundaryConditionDiscretizationSpectralConform::
00723 transposedTimesX(const BaseVector& X, BaseVector& Z) const
00724 {
00725 const Vector<real_t>& u = dynamic_cast<const Vector<real_t>&>(X);
00726 Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(Z);
00727
00728 switch(this->problem().type()) {
00729 case Problem::variational: {
00730 const VariationalProblem& variationalProblem
00731 = dynamic_cast<const VariationalProblem&>(this->problem());
00732
00733 for (VariationalProblem::bilinearBorderOperatorConst_iterator
00734 iBorderOperator = variationalProblem.beginBilinearBorderOperator();
00735 iBorderOperator != variationalProblem.endBilinearBorderOperator();
00736 ++iBorderOperator) {
00737
00738 const VariationalBilinearBorderOperator& borderOperator = **iBorderOperator;
00739 switch (borderOperator.type()) {
00740 case VariationalBilinearBorderOperator::alphaUV: {
00741
00742 const VariationalBorderOperatorAlphaUV& operatorAlphaUV
00743 = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(borderOperator);
00744 ConstReferenceCounting<Boundary> boundary = operatorAlphaUV.boundary();
00745
00746 switch(boundary->type()) {
00747 case Boundary::references: {
00748 const BoundaryReferences& boundaryReference = dynamic_cast<const BoundaryReferences&>(*boundary);
00749 for (BoundaryReferences::ReferencesSet::const_iterator
00750 iBoundary = boundaryReference.references().begin();
00751 iBoundary != boundaryReference.references().end(); ++iBoundary) {
00752
00753 const ScalarFunctionBase& alpha = operatorAlphaUV.alpha();
00754 const size_t boundaryNumber = *iBoundary;
00755
00756 const size_t& numberOfUnknown = this->problem().numberOfUnknown();
00757
00758 Vector<Vector<Vector<real_t> > > u_jk(numberOfUnknown);
00759 Vector<Vector<Vector<real_t> > > u_js(numberOfUnknown);
00760 Vector<Vector<Vector<real_t> > > u_rs(numberOfUnknown);
00761 Vector<Vector<Vector<real_t> > > v_rp(numberOfUnknown);
00762 Vector<Vector<Vector<real_t> > > v_np(numberOfUnknown);
00763
00764 TinyVector<3, size_t> direction;
00765 for (size_t i=0; i<3; ++i) {
00766 direction[i] = (boundaryNumber/2+i+1) % 3;
00767 }
00768
00769 Vector<TinyVector<3,Vector<Vector<real_t> > > > quadratureBaseValuesU(numberOfUnknown);
00770 for (size_t unknownNumber =0; unknownNumber < numberOfUnknown; ++unknownNumber){
00771
00772 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >&unknownBasis
00773 = __basisU[unknownNumber];
00774
00775 for (size_t i=0; i<2; ++i) {
00776 quadratureBaseValuesU[unknownNumber][i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00777
00778 for (size_t j=0; j < (*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00779
00780 real_t xj=(*__transformU[unknownNumber][direction[i]]).inverse((*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j)));
00781
00782 quadratureBaseValuesU[unknownNumber][i][j].resize(unknownBasis[direction[i]]->dimension());
00783 unknownBasis[direction[i]]->getValues(xj,quadratureBaseValuesU[unknownNumber][i][j]);
00784 }
00785 }
00786
00787 quadratureBaseValuesU[unknownNumber][2].resize(2);
00788 for (size_t j=0; j<2 ; ++j) {
00789 real_t xj= -1. + 2.*j;
00790 quadratureBaseValuesU[unknownNumber][2][j].resize(unknownBasis[direction[2]]->dimension());
00791 unknownBasis[direction[2]]->getValues(xj,quadratureBaseValuesU[unknownNumber][2][j]);
00792 }
00793 }
00794
00795
00796
00797 TinyVector<2, Vector<real_t> >nodes;
00798 for (size_t i=0; i<2; ++i) {
00799 nodes[i].resize((*__gaussLobatto[direction[i]]).numberOfPoints());
00800 for (size_t j=0; j<(*__gaussLobatto[direction[i]]).numberOfPoints(); ++j) {
00801 nodes[i][j] =(*__transform[direction[i]])((*__gaussLobatto[direction[i]])(j));
00802 }
00803 }
00804 const real_t jacobianDet = __transform[direction[0]]->determinant()*
00805 __transform[direction[1]]->determinant();
00806
00807
00808 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >&unknownBasis
00809 = __basisU[operatorAlphaUV.unknownNumber()];
00810
00811 u_jk[operatorAlphaUV.unknownNumber()].resize(unknownBasis[direction[0]]->dimension());
00812 for(size_t j=0; j< unknownBasis[direction[0]]->dimension(); ++j){
00813 u_jk[operatorAlphaUV.unknownNumber()][j].resize( unknownBasis[direction[1]]->dimension());
00814 for(size_t k=0; k < unknownBasis[direction[1]]->dimension(); ++k){
00815 u_jk[operatorAlphaUV.unknownNumber()][j][k]=0;
00816 }
00817 }
00818
00819 {
00820 const TinyVector<3,size_t>& dofShape = __dofShape[operatorAlphaUV.unknownNumber()];
00821 TinyVector<3,size_t> J;
00822 for(size_t j=0; j< unknownBasis[direction[0]]->dimension(); ++j){
00823 J[direction[0]]=j;
00824 for(size_t k=0; k < unknownBasis[direction[1]]->dimension(); ++k){
00825 J[direction[1]]=k;
00826 real_t& u_jk_value = u_jk[operatorAlphaUV.unknownNumber()][j][k];
00827 for(size_t i=0; i < unknownBasis[direction[2]]->dimension(); ++i){
00828 J[direction[2]]=i;
00829 const size_t dofNumberU =
00830 J[0] *dofShape[1]*dofShape[2] + J[1] *dofShape[2] +J[2] ;
00831 u_jk_value
00832 += quadratureBaseValuesU[operatorAlphaUV.unknownNumber()][2][boundaryNumber%2][i]
00833 *u[__degreeOfFreedomSet(operatorAlphaUV.testFunctionNumber(),dofNumberU)];
00834 }
00835 }
00836 }
00837 }
00838
00839 u_js[operatorAlphaUV.unknownNumber()].resize(unknownBasis[direction[0]]->dimension());
00840 for(size_t j=0; j< unknownBasis[direction[0]]->dimension(); ++j){
00841 u_js[operatorAlphaUV.unknownNumber()][j].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00842 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00843 u_js[operatorAlphaUV.unknownNumber()][j][s]=0;
00844 }
00845 }
00846
00847 for(size_t j=0; j < unknownBasis[direction[0]]->dimension(); ++j){
00848 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00849 real_t& u_js_value = u_js[operatorAlphaUV.unknownNumber()][j][s];
00850 for(size_t k=0; k < unknownBasis[direction[1]]->dimension(); ++k){
00851 u_js_value
00852 += quadratureBaseValuesU[operatorAlphaUV.unknownNumber()][1][s][k]
00853 * u_jk[operatorAlphaUV.unknownNumber()][j][k];
00854 }
00855 }
00856 }
00857 u_rs[operatorAlphaUV.unknownNumber()].resize(__gaussLobatto[direction[0]]->numberOfPoints());
00858 for(size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00859 u_rs[operatorAlphaUV.unknownNumber()][r].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00860 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00861 u_rs[operatorAlphaUV.unknownNumber()][r][s]=0;
00862 }
00863 }
00864
00865 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00866 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00867 real_t& u_rs_value = u_rs[operatorAlphaUV.unknownNumber()][r][s];
00868 for(size_t j=0; j < unknownBasis[direction[0]]->dimension(); ++j){
00869 u_rs_value
00870 += quadratureBaseValuesU[operatorAlphaUV.unknownNumber()][0][r][j]
00871 *u_js[operatorAlphaUV.unknownNumber()][j][s];
00872 }
00873 }
00874 }
00875
00876 Vector<Vector<real_t> > alpha_rs(__gaussLobatto[direction[0]]->numberOfPoints());
00877 for(size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00878 alpha_rs[r].resize(__gaussLobatto[direction[1]]->numberOfPoints());
00879 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00880 alpha_rs[r][s]=0;
00881 }
00882 }
00883
00884 for (size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r) {
00885 const real_t& x = nodes[0][r];
00886 for (size_t s=0; s <__gaussLobatto[direction[1]]->numberOfPoints(); ++s) {
00887 const real_t& y = nodes[1][s];
00888
00889 switch(boundaryNumber) {
00890 case 0: {
00891 alpha_rs[r][s] = alpha(__mesh.shape().a(0), x, y);
00892 break;
00893 }
00894 case 1: {
00895 alpha_rs[r][s] = alpha(__mesh.shape().b(0), x, y);
00896 break;
00897 }
00898 case 2: {
00899 alpha_rs[r][s] = alpha(y, __mesh.shape().a(1), x);
00900 break;
00901 }
00902 case 3: {
00903 alpha_rs[r][s] = alpha(y, __mesh.shape().b(1), x);
00904 break;
00905 }
00906 case 4: {
00907 alpha_rs[r][s] = alpha(x, y, __mesh.shape().a(2));
00908 break;
00909 }
00910 case 5: {
00911 alpha_rs[r][s] = alpha(x, y, __mesh.shape().b(2));
00912 break;
00913 }
00914 }
00915 }
00916 }
00917
00918 const TinyVector<3, ConstReferenceCounting<LegendreBasis> >&testBasis
00919 = __basisU[operatorAlphaUV.testFunctionNumber()];
00920
00921 v_rp[operatorAlphaUV.testFunctionNumber()].resize(__gaussLobatto[direction[0]]->numberOfPoints());
00922
00923 for(size_t r=0; r<__gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00924 v_rp[operatorAlphaUV.testFunctionNumber()][r].resize(testBasis[direction[1]]->dimension());
00925 for(size_t p=0; p < testBasis[direction[1]]->dimension(); ++p){
00926 v_rp[operatorAlphaUV.testFunctionNumber()][r][p]=0;
00927 }
00928 }
00929
00930 for(size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00931 for(size_t p=0; p < testBasis[direction[1]]->dimension(); ++p){
00932 real_t& v_rp_value = v_rp[operatorAlphaUV.testFunctionNumber()][r][p];
00933 for(size_t s=0; s<__gaussLobatto[direction[1]]->numberOfPoints(); ++s){
00934 v_rp_value
00935 += quadratureBaseValuesU[operatorAlphaUV.testFunctionNumber()][1][s][p]
00936 *u_rs[operatorAlphaUV.unknownNumber()][r][s]
00937 *__gaussLobatto[direction[1]]->weight(s)
00938 * alpha_rs[r][s]*jacobianDet;
00939 }
00940 }
00941 }
00942
00943 v_np[operatorAlphaUV.testFunctionNumber()].resize(testBasis[direction[0]]->dimension());
00944 for(size_t n=0; n < testBasis[direction[0]]->dimension(); ++n){
00945 v_np[operatorAlphaUV.testFunctionNumber()][n].resize(testBasis[direction[1]]->dimension());
00946 for(size_t p=0; p < testBasis[direction[1]]->dimension(); ++p){
00947 v_np[operatorAlphaUV.testFunctionNumber()][n][p]=0;
00948 }
00949 }
00950
00951 for(size_t n=0; n < testBasis[direction[0]]->dimension(); ++n){
00952 for(size_t p=0; p < testBasis[direction[1]]->dimension(); ++p){
00953 real_t& v_np_value = v_np[operatorAlphaUV.testFunctionNumber()][n][p];
00954 for(size_t r=0; r < __gaussLobatto[direction[0]]->numberOfPoints(); ++r){
00955 v_np_value
00956 += quadratureBaseValuesU[operatorAlphaUV.testFunctionNumber()][0][r][n]
00957 *v_rp[operatorAlphaUV.testFunctionNumber()][r][p]*
00958 __gaussLobatto[direction[0]]->weight(r);
00959 }
00960 }
00961 }
00962
00963 const TinyVector<3,size_t>& dofShapeV = __dofShape[operatorAlphaUV.testFunctionNumber()];
00964 TinyVector<3,size_t> I;
00965
00966 for(size_t n=0; n < testBasis[direction[0]]->dimension(); ++n){
00967 I[direction[0]]=n;
00968 for(size_t p=0; p < testBasis[direction[1]]->dimension(); ++p){
00969 I[direction[1]]=p;
00970 for(size_t m=0; m < testBasis[direction[2]]->dimension(); ++m){
00971 I[direction[2]]=m;
00972 const size_t dofNumberV =
00973 I[0]* dofShapeV[1]*dofShapeV[2] + I[1]*dofShapeV[2] + I[2];
00974
00975 v[__degreeOfFreedomSet(operatorAlphaUV.unknownNumber(),dofNumberV)]
00976 += quadratureBaseValuesU[operatorAlphaUV.testFunctionNumber()][2][boundaryNumber%2][m]
00977 *v_np[operatorAlphaUV.testFunctionNumber()][n][p];
00978 }
00979 }
00980 }
00981 }
00982 break;
00983 }
00984 default: {
00985 throw ErrorHandler(__FILE__,__LINE__,
00986 "not implemented",
00987 ErrorHandler::unexpected);
00988 }
00989 }
00990 break;
00991 }
00992 default: {
00993 throw ErrorHandler(__FILE__,__LINE__,
00994 "not implemented",
00995 ErrorHandler::unexpected);
00996 }
00997 }
00998 }
00999 break;
01000 }
01001 case Problem::pde: {
01002 break;
01003 }
01004 default: {
01005 throw ErrorHandler(__FILE__,__LINE__,
01006 "not implemented",
01007 ErrorHandler::unexpected);
01008 }
01009 }
01010 }
01011
01012
01013
01014
01015
01016
