BoundaryConditionDiscretizationSpectralConform Class Reference

#include <BoundaryConditionDiscretizationSpectralConform.hpp>

Inheritance diagram for BoundaryConditionDiscretizationSpectralConform:

Inheritance graph
[legend]
Collaboration diagram for BoundaryConditionDiscretizationSpectralConform:

Collaboration graph
[legend]

List of all members.

Public Member Functions

void getDiagonal (BaseVector &X) const
void setMatrix (ReferenceCounting< BaseMatrix > A, ReferenceCounting< BaseVector > b) const
void setSecondMember (ReferenceCounting< BaseMatrix > A, ReferenceCounting< BaseVector > b) const
void timesX (const BaseVector &X, BaseVector &Z) const
void transposedTimesX (const BaseVector &X, BaseVector &Z) const
 BoundaryConditionDiscretizationSpectralConform (const Problem &problem, const SpectralMesh &mesh, const DegreeOfFreedomSet &dof, const DiscretizationType &discretizationType)
 ~BoundaryConditionDiscretizationSpectralConform ()
const bool & dirichlet (const size_t i) const
const real_t & dirichletValue (const size_t i) const
const Vector< bool > & getDirichletList () const
const Problemproblem () const

Protected Attributes

const Problem__problem
Vector< bool > __dirichletList
 true for vertices supporting Dirichlet.
Vector< real_t > __dirichletValues

Private Attributes

const DegreeOfFreedomSet__degreeOfFreedomSet
const SpectralMesh__mesh
Vector< TinyVector< 3, size_t > > __dofShape
TinyVector
< 3, ConstReferenceCounting
< Interval > > 
__interval
Vector< TinyVector
< 3, ConstReferenceCounting
< Interval > > > 
__intervalU
TinyVector
< 3, ConstReferenceCounting
< SpectralConformTransformation > > 
__transform
Vector< TinyVector
< 3, ConstReferenceCounting
< SpectralConformTransformation > > > 
__transformU
TinyVector
< 3, ConstReferenceCounting
< GaussLobatto > > 
__gaussLobatto
Vector< TinyVector
< 3, ConstReferenceCounting
< LegendreBasis > > > 
__basisU


Detailed Description

Definition at line 44 of file BoundaryConditionDiscretizationSpectralConform.hpp.


Constructor & Destructor Documentation

BoundaryConditionDiscretizationSpectralConform::BoundaryConditionDiscretizationSpectralConform ( const Problem problem,
const SpectralMesh mesh,
const DegreeOfFreedomSet dof,
const DiscretizationType discretizationType 
) [inline]

Definition at line 75 of file BoundaryConditionDiscretizationSpectralConform.hpp.

References __basisU, __dofShape, __intervalU, __transformU, ScalarDiscretizationTypeSpectral::a(), ScalarDiscretizationTypeSpectral::b(), ScalarDiscretizationTypeSpectral::degrees(), ScalarDiscretizationTypeBase::name(), Problem::numberOfUnknown(), BoundaryConditionDiscretization::problem(), ScalarDiscretizationTypeBase::spectralLegendre, and ErrorHandler::unexpected.

00079     : BoundaryConditionDiscretization(problem,dof),
00080       __degreeOfFreedomSet(dof),
00081       __mesh(mesh),
00082       __dofShape(this->problem().numberOfUnknown()),
00083       __interval(new Interval(__mesh.shape().a()[0],__mesh.shape().b()[0]),
00084                  new Interval(__mesh.shape().a()[1],__mesh.shape().b()[1]),
00085                  new Interval(__mesh.shape().a()[2],__mesh.shape().b()[2])),
00086       
00087       __intervalU(this->problem().numberOfUnknown()),
00088       __transform(new SpectralConformTransformation(*__interval[0]),
00089                   new SpectralConformTransformation(*__interval[1]),
00090                   new SpectralConformTransformation(*__interval[2])),
00091       
00092       __transformU(this->problem().numberOfUnknown()),
00093       
00094       __gaussLobatto(GaussLobattoManager::instance().getReference(__mesh.degree(0)+1),
00095                      GaussLobattoManager::instance().getReference(__mesh.degree(1)+1),
00096                      GaussLobattoManager::instance().getReference(__mesh.degree(2)+1)),
00097       __basisU(this->problem().numberOfUnknown())
00098   {
00099     for (size_t i=0; i < this->problem().numberOfUnknown(); i++) {
00100       if (discretizationType[i].type() != ScalarDiscretizationTypeBase::spectralLegendre) {
00101         throw ErrorHandler(__FILE__,__LINE__,
00102                            "discretization '"
00103                            +ScalarDiscretizationTypeBase::name(discretizationType[i])
00104                            +"' is incompatible with spectral method",
00105                            ErrorHandler::unexpected);
00106       }
00107       const ScalarDiscretizationTypeSpectral& discretization
00108         = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(discretizationType[i]);
00109       
00110       __dofShape[i] = discretization.degrees()+TinyVector<3,size_t>(1,1,1);
00111       
00112       __intervalU[i][0]= new Interval(discretization.a()[0],discretization.b()[0]);
00113       __intervalU[i][1]= new Interval(discretization.a()[1],discretization.b()[1]);
00114       __intervalU[i][2]= new Interval(discretization.a()[2],discretization.b()[2]);
00115       __transformU[i][0] = new SpectralConformTransformation(*__intervalU[i][0]);
00116       __transformU[i][1] = new SpectralConformTransformation(*__intervalU[i][1]);
00117       __transformU[i][2] = new SpectralConformTransformation(*__intervalU[i][2]); 
00118       __basisU[i][0]= new LegendreBasis(discretization.degrees()[0]);
00119       __basisU[i][1]= new LegendreBasis(discretization.degrees()[1]);
00120       __basisU[i][2]= new LegendreBasis(discretization.degrees()[2]);
00121     }
00122     
00123   }

Here is the call graph for this function:

BoundaryConditionDiscretizationSpectralConform::~BoundaryConditionDiscretizationSpectralConform (  )  [inline]

Definition at line 125 of file BoundaryConditionDiscretizationSpectralConform.hpp.

00126   {
00127     ;
00128   }


Member Function Documentation

void BoundaryConditionDiscretizationSpectralConform::getDiagonal ( BaseVector X  )  const [virtual]

Implements BoundaryConditionDiscretization.

Definition at line 43 of file BoundaryConditionDiscretizationSpectralConform.cpp.

References __basisU, __degreeOfFreedomSet, __dofShape, __gaussLobatto, __mesh, __transform, __transformU, Structured3DMeshShape::a(), VariationalBilinearBorderOperator::alphaUV, Structured3DMeshShape::b(), VariationalProblem::beginBilinearBorderOperator(), VariationalProblem::endBilinearBorderOperator(), Problem::numberOfUnknown(), Problem::pde, BoundaryConditionDiscretization::problem(), BoundaryReferences::references(), Boundary::references, Vector< T >::resize(), SpectralMesh::shape(), VariationalBorderOperator::testFunctionNumber(), VariationalBilinearBorderOperator::type(), Problem::type(), ErrorHandler::unexpected, VariationalBilinearBorderOperator::unknownNumber(), and Problem::variational.

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 }

Here is the call graph for this function:

void BoundaryConditionDiscretizationSpectralConform::setMatrix ( ReferenceCounting< BaseMatrix A,
ReferenceCounting< BaseVector b 
) const [virtual]

Implements BoundaryConditionDiscretization.

Definition at line 231 of file BoundaryConditionDiscretizationSpectralConform.cpp.

References BaseMatrix::doubleHashedMatrix, UnAssembledMatrix::setBoundaryConditions(), BaseMatrix::unAssembled, and ErrorHandler::unexpected.

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 }

Here is the call graph for this function:

void BoundaryConditionDiscretizationSpectralConform::setSecondMember ( ReferenceCounting< BaseMatrix A,
ReferenceCounting< BaseVector b 
) const [virtual]

Implements BoundaryConditionDiscretization.

Definition at line 257 of file BoundaryConditionDiscretizationSpectralConform.cpp.

References __basisU, __degreeOfFreedomSet, __dofShape, __gaussLobatto, __mesh, __transform, __transformU, Structured3DMeshShape::a(), Structured3DMeshShape::b(), VariationalProblem::beginLinearBorderOperator(), VariationalProblem::endLinearBorderOperator(), ffout(), VariationalLinearBorderOperator::FV, Problem::numberOfUnknown(), BoundaryConditionDiscretization::problem(), BoundaryReferences::references(), Boundary::references, Vector< T >::resize(), SpectralMesh::shape(), VariationalBorderOperator::testFunctionNumber(), VariationalLinearBorderOperator::type(), Problem::type(), ErrorHandler::unexpected, and Problem::variational.

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             // nodes construction           
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 }

Here is the call graph for this function:

void BoundaryConditionDiscretizationSpectralConform::timesX ( const BaseVector X,
BaseVector Z 
) const [virtual]

Implements BoundaryConditionDiscretization.

Definition at line 439 of file BoundaryConditionDiscretizationSpectralConform.cpp.

References __basisU, __degreeOfFreedomSet, __dofShape, __gaussLobatto, __mesh, __transform, __transformU, Structured3DMeshShape::a(), VariationalBilinearBorderOperator::alphaUV, Structured3DMeshShape::b(), VariationalProblem::beginBilinearBorderOperator(), VariationalProblem::endBilinearBorderOperator(), Problem::numberOfUnknown(), Problem::pde, BoundaryConditionDiscretization::problem(), BoundaryReferences::references(), Boundary::references, Vector< T >::resize(), SpectralMesh::shape(), VariationalBorderOperator::testFunctionNumber(), VariationalBilinearBorderOperator::type(), Problem::type(), ErrorHandler::unexpected, VariationalBilinearBorderOperator::unknownNumber(), and Problem::variational.

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             // nodes construction
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 }

Here is the call graph for this function:

void BoundaryConditionDiscretizationSpectralConform::transposedTimesX ( const BaseVector X,
BaseVector Z 
) const [virtual]

Implements BoundaryConditionDiscretization.

Definition at line 723 of file BoundaryConditionDiscretizationSpectralConform.cpp.

References __basisU, __degreeOfFreedomSet, __dofShape, __gaussLobatto, __mesh, __transform, __transformU, Structured3DMeshShape::a(), VariationalBilinearBorderOperator::alphaUV, Structured3DMeshShape::b(), VariationalProblem::beginBilinearBorderOperator(), VariationalProblem::endBilinearBorderOperator(), Problem::numberOfUnknown(), Problem::pde, BoundaryConditionDiscretization::problem(), BoundaryReferences::references(), Boundary::references, Vector< T >::resize(), SpectralMesh::shape(), VariationalBorderOperator::testFunctionNumber(), VariationalBilinearBorderOperator::type(), Problem::type(), ErrorHandler::unexpected, VariationalBilinearBorderOperator::unknownNumber(), and Problem::variational.

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             //construction des nodes
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             //const size_t unknownNumber= operatorAlphaUV.unknownNumber();
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 }

Here is the call graph for this function:

const bool& BoundaryConditionDiscretization::dirichlet ( const size_t  i  )  const [inline, inherited]

Definition at line 54 of file BoundaryConditionDiscretization.hpp.

References BoundaryConditionDiscretization::__dirichletList.

00055   {
00056     return __dirichletList[i];
00057   }

const real_t& BoundaryConditionDiscretization::dirichletValue ( const size_t  i  )  const [inline, inherited]

Definition at line 59 of file BoundaryConditionDiscretization.hpp.

References BoundaryConditionDiscretization::__dirichletValues.

00060   {
00061     return __dirichletValues[i];
00062   }

const Vector<bool>& BoundaryConditionDiscretization::getDirichletList (  )  const [inline, inherited]

Definition at line 64 of file BoundaryConditionDiscretization.hpp.

References BoundaryConditionDiscretization::__dirichletList.

00065   {
00066     return __dirichletList;
00067   }

const Problem& BoundaryConditionDiscretization::problem (  )  const [inline, inherited]


Member Data Documentation

const Problem& BoundaryConditionDiscretization::__problem [protected, inherited]

Vector<bool> BoundaryConditionDiscretization::__dirichletList [mutable, protected, inherited]

Vector<real_t> BoundaryConditionDiscretization::__dirichletValues [mutable, protected, inherited]


The documentation for this class was generated from the following files:

Generated on Wed Nov 19 00:04:35 2008 for FreeFEM3D (aka ff3d) by  doxygen 1.5.6