00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef SPECTRAL_FUNCTION_HPP
00022 #define SPECTRAL_FUNCTION_HPP
00023
00024 #include <ScalarFunctionBase.hpp>
00025
00026
00027 #include <SpectralMesh.hpp>
00028 #include <Vector.hpp>
00029 #include <Interval.hpp>
00030 #include <LegendreBasis.hpp>
00031
00032 #include <GaussLobattoManager.hpp>
00033
00034 #include <SpectralConformTransformation.hpp>
00035
00036 #include <ScalarFunctionConstant.hpp>
00037
00046 class SpectralFunction
00047 : public ScalarFunctionBase
00048 {
00049 private:
00050 Vector<real_t> __values;
00053 ConstReferenceCounting<SpectralMesh>
00054 __mesh;
00057 const ScalarDiscretizationTypeBase::Type
00058 __discretizationType;
00060 const Interval __intervalX;
00061 const Interval __intervalY;
00062 const Interval __intervalZ;
00064 const SpectralConformTransformation
00065 __transformX;
00067 const SpectralConformTransformation
00068 __transformY;
00070 const SpectralConformTransformation
00071 __transformZ;
00073 const GaussLobatto&
00074 __gaussLobattoX;
00075 const GaussLobatto&
00076 __gaussLobattoY;
00077 const GaussLobatto&
00078 __gaussLobattoZ;
00080 const LegendreBasis __basisX;
00081 const LegendreBasis __basisY;
00082 const LegendreBasis __basisZ;
00084 TinyVector<3, ConstReferenceCounting<Interval> > __interval;
00085 TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > __transform;
00086 TinyVector<3, ConstReferenceCounting<GaussLobatto> > __gaussLobatto;
00087 TinyVector<3, ConstReferenceCounting<LegendreBasis> > __basis;
00088
00096 std::ostream& __put(std::ostream& os) const
00097 {
00098 os << "{sfunction}";
00099 return os;
00100 }
00101
00102 public:
00108 const Vector<real_t>& values() const
00109 {
00110 return __values;
00111 }
00112
00120 inline real_t& operator[](const size_t& i)
00121 {
00122 return __values[i];
00123 }
00124
00132 inline const real_t& operator[](const size_t& i) const
00133 {
00134 return __values[i];
00135 }
00136
00142 ConstReferenceCounting<SpectralMesh> mesh() const
00143 {
00144 return __mesh;
00145 }
00146
00152 bool canBeSimplified() const
00153 {
00154 return false;
00155 }
00156
00162 const ScalarDiscretizationTypeBase::Type& discretizationType() const
00163 {
00164 return __discretizationType;
00165 }
00166
00174 real_t operator()(const TinyVector<3,real_t>& x) const
00175 {
00176
00177 Vector<real_t> baseValuesX(__basisX.dimension());
00178 Vector<real_t> baseValuesY(__basisY.dimension());
00179 Vector<real_t> baseValuesZ(__basisZ.dimension());
00180
00181 __basisX.getValues(__transformX.inverse(x[0]),baseValuesX);
00182 __basisY.getValues(__transformY.inverse(x[1]),baseValuesY);
00183 __basisZ.getValues(__transformZ.inverse(x[2]),baseValuesZ);
00184
00185 size_t l = 0;
00186 real_t sum = 0;
00187 for (size_t i=0; i<__basisX.dimension(); ++i) {
00188 const real_t vi = baseValuesX[i];
00189 for (size_t j=0;j<__basisY.dimension(); ++j) {
00190 const real_t vj = baseValuesY[j];
00191 const real_t vi_vj = vi*vj;
00192 for (size_t k=0; k<__basisZ.dimension(); ++k) {
00193 const real_t vk = baseValuesZ[k];
00194 const real_t vi_vj_vk = vi_vj * vk;
00195 sum += vi_vj_vk*__values[l];
00196 l++;
00197 }
00198 }
00199 }
00200
00201 return sum;
00202 }
00203
00209 void operator=(const ScalarFunctionBase& f)
00210 {
00211 __values = 0;
00212
00217 Vector<real_t> nodesX(__gaussLobattoX.numberOfPoints());
00218 for (size_t i=0; i<__gaussLobattoX.numberOfPoints(); ++i) {
00219 nodesX[i] = __transformX(__gaussLobattoX(i));
00220 }
00221
00222 Vector<real_t> nodesY(__gaussLobattoY.numberOfPoints());
00223 for (size_t i=0; i<__gaussLobattoY.numberOfPoints(); ++i) {
00224 nodesY[i] = __transformY(__gaussLobattoY(i));
00225 }
00226
00227 Vector<real_t> nodesZ(__gaussLobattoZ.numberOfPoints());
00228 for (size_t i=0; i<__gaussLobattoZ.numberOfPoints(); ++i) {
00229 nodesZ[i] = __transformZ(__gaussLobattoZ(i));
00230 }
00231
00236 Vector<Vector<Vector<real_t> > > f_rst(__gaussLobattoX.numberOfPoints());
00237 for (size_t r=0; r<__gaussLobattoX.numberOfPoints(); ++r) {
00238 f_rst[r].resize(__gaussLobattoY.numberOfPoints());
00239 for (size_t s=0; s<__gaussLobattoY.numberOfPoints(); ++s) {
00240 f_rst[r][s].resize(__gaussLobattoZ.numberOfPoints());
00241 f_rst[r][s] = 0;
00242 }
00243 }
00244
00245 for (size_t r=0; r < __gaussLobattoX.numberOfPoints(); ++r) {
00246 const real_t& x = nodesX[r];
00247 for (size_t s=0; s < __gaussLobattoY.numberOfPoints(); ++s) {
00248 const real_t& y = nodesY[s];
00249 for (size_t t=0; t < __gaussLobattoZ.numberOfPoints(); ++t) {
00250 const real_t& z = nodesZ[t];
00251 f_rst[r][s][t] = f(x,y,z);
00252 }
00253 }
00254 }
00255
00260 Vector<Vector<real_t> > quadratureBaseValuesX(__gaussLobattoX.numberOfPoints());
00261 for (size_t i=0; i < __gaussLobattoX.numberOfPoints(); ++i) {
00262 real_t xi= __gaussLobattoX(i);
00263 quadratureBaseValuesX[i].resize(__basisX.dimension());
00264 __basisX.getValues(xi,quadratureBaseValuesX[i]);
00265 }
00266
00267 Vector<Vector<real_t> > quadratureBaseValuesY(__gaussLobattoY.numberOfPoints());
00268 for (size_t j=0; j < __gaussLobattoY.numberOfPoints(); ++j) {
00269 real_t yj= __gaussLobattoY(j);
00270 quadratureBaseValuesY[j].resize(__basisY.dimension());
00271 __basisY.getValues(yj,quadratureBaseValuesY[j]);
00272 }
00273
00274 Vector<Vector<real_t> > quadratureBaseValuesZ(__gaussLobattoZ.numberOfPoints());
00275 for (size_t k=0; k < __gaussLobattoZ.numberOfPoints(); ++k) {
00276 real_t zk= __gaussLobattoZ(k);
00277 quadratureBaseValuesZ[k].resize(__basisZ.dimension());
00278 __basisZ.getValues(zk,quadratureBaseValuesZ[k]);
00279 }
00280
00285 Vector<Vector<Vector<real_t> > > values_rsk(__gaussLobattoX.numberOfPoints());
00286 for (size_t r=0; r<__gaussLobattoX.numberOfPoints(); ++r) {
00287 values_rsk[r].resize( __gaussLobattoY.numberOfPoints());
00288 for (size_t s=0; s<__gaussLobattoY.numberOfPoints(); ++s) {
00289 values_rsk[r][s].resize( __basisZ.dimension());
00290 values_rsk[r][s] = 0;
00291 }
00292 }
00293
00294 {
00295 for (size_t r=0; r < __gaussLobattoX.numberOfPoints(); ++r) {
00296 for (size_t s=0; s < __gaussLobattoY.numberOfPoints(); ++s) {
00297 for (size_t t=0; t < __gaussLobattoZ.numberOfPoints(); ++t) {
00298 const real_t& wt = __gaussLobattoZ.weight(t);
00299 const real_t value_rst = f_rst[r][s][t]* wt;
00300 for (size_t k=0; k<__basisZ.dimension(); ++k) {
00301 values_rsk[r][s][k] += value_rst * quadratureBaseValuesZ[t][k];
00302 }
00303 }
00304 }
00305 }
00306 }
00307
00308 Vector<Vector<Vector<real_t> > > values_rjk(__gaussLobattoX.numberOfPoints());
00309 for (size_t r=0; r<__gaussLobattoX.numberOfPoints(); ++r) {
00310 values_rjk[r].resize( __basisY.dimension());
00311 for (size_t j=0; j<__basisY.dimension(); ++j) {
00312 values_rjk[r][j].resize( __basisZ.dimension());
00313 values_rjk[r][j] = 0;
00314 }
00315 }
00316
00317 {
00318 for (size_t r=0; r < __gaussLobattoX.numberOfPoints(); ++r) {
00319 for (size_t s=0; s < __gaussLobattoY.numberOfPoints(); ++s) {
00320 const real_t& ws = __gaussLobattoY.weight(s);
00321 for (size_t k=0; k < __basisZ.dimension(); ++k) {
00322 const real_t value_rsk = values_rsk[r][s][k]* ws;
00323 for (size_t j=0; j<__basisY.dimension(); ++j) {
00324 values_rjk[r][j][k] += value_rsk * quadratureBaseValuesY[s][j];
00325 }
00326 }
00327 }
00328 }
00329 }
00330
00331 {
00332
00333 for (size_t r=0; r < __gaussLobattoX.numberOfPoints(); ++r) {
00334 const real_t& wr = __gaussLobattoX.weight(r);
00335
00336 for (size_t j=0; j < __basisY.dimension(); ++j) {
00337 for (size_t k=0; k < __basisZ.dimension(); ++k) {
00338 const real_t value_rjk = values_rjk[r][j][k]* wr;
00339 for (size_t i=0; i<__basisX.dimension(); ++i) {
00340 __values[__mesh->dofNumber(i,j,k)] +=
00341 value_rjk * quadratureBaseValuesX[r][i];
00342
00343 }
00344 }
00345 }
00346 }
00347 }
00348
00352
00353 for (size_t i=0; i<__basisX.dimension(); ++i) {
00354 const real_t wi_8 = 1./8. * (2*i+1);
00355 for (size_t j=0;j<__basisY.dimension(); ++j) {
00356 const real_t wi_wj_8 = wi_8 * (2*j+1);
00357 for (size_t k=0; k<__basisZ.dimension(); ++k) {
00358 const real_t wi_wj_wk_8 = wi_wj_8 * (2*k+1);
00359 __values[__mesh->dofNumber(i,j,k)] *= wi_wj_wk_8;
00360
00361 }
00362 }
00363 }
00364 }
00365
00373 real_t dx(const TinyVector<3>& x) const
00374 {
00375 Vector<real_t> baseDerivativeValuesX(__basisX.dimension());
00376 Vector<real_t> baseValuesY(__basisY.dimension());
00377 Vector<real_t> baseValuesZ(__basisZ.dimension());
00378
00379 __basisX.getDerivativeValues(__transformX.inverse(x[0]),baseDerivativeValuesX);
00380 __basisY.getValues(__transformY.inverse(x[1]),baseValuesY);
00381 __basisZ.getValues(__transformZ.inverse(x[2]),baseValuesZ);
00382
00383 size_t l = 0;
00384 real_t sum = 0;
00385 for (size_t i=0; i<__basisX.dimension(); ++i) {
00386 const real_t vi = baseDerivativeValuesX[i];
00387 for (size_t j=0;j<__basisY.dimension(); ++j) {
00388 const real_t vj = baseValuesY[j];
00389 const real_t vi_vj = vi*vj;
00390 for (size_t k=0; k<__basisZ.dimension(); ++k) {
00391 const real_t vk = baseValuesZ[k];
00392 const real_t vi_vj_vk = vi_vj * vk;
00393 sum += vi_vj_vk*__values[l];
00394 l++;
00395 }
00396 }
00397 }
00398
00399 return sum * __transformX.inverseDeterminant();
00400 }
00401
00409 real_t dy(const TinyVector<3>& x) const
00410 {
00411 Vector<real_t> baseValuesX(__basisX.dimension());
00412 Vector<real_t> baseDerivativeValuesY(__basisY.dimension());
00413 Vector<real_t> baseValuesZ(__basisZ.dimension());
00414
00415 __basisX.getValues(__transformX.inverse(x[0]),baseValuesX);
00416 __basisY.getDerivativeValues(__transformY.inverse(x[1]),baseDerivativeValuesY);
00417 __basisZ.getValues(__transformZ.inverse(x[2]),baseValuesZ);
00418
00419 size_t l = 0;
00420 real_t sum = 0;
00421 for (size_t i=0; i<__basisX.dimension(); ++i) {
00422 const real_t vi = baseValuesX[i];
00423 for (size_t j=0;j<__basisY.dimension(); ++j) {
00424 const real_t vj = baseDerivativeValuesY[j];
00425 const real_t vi_vj = vi*vj;
00426 for (size_t k=0; k<__basisZ.dimension(); ++k) {
00427 const real_t vk = baseValuesZ[k];
00428 const real_t vi_vj_vk = vi_vj * vk;
00429 sum += vi_vj_vk*__values[l];
00430 l++;
00431 }
00432 }
00433 }
00434
00435 return sum * __transformY.inverseDeterminant();
00436 }
00437
00445 real_t dz(const TinyVector<3>& x) const
00446 {
00447 Vector<real_t> baseValuesX(__basisX.dimension());
00448 Vector<real_t> baseValuesY(__basisY.dimension());
00449 Vector<real_t> baseDerivativeValuesZ(__basisZ.dimension());
00450
00451 __basisX.getValues(__transformX.inverse(x[0]),baseValuesX);
00452 __basisY.getValues(__transformY.inverse(x[1]),baseValuesY);
00453 __basisZ.getDerivativeValues(__transformZ.inverse(x[2]),baseDerivativeValuesZ);
00454
00455 size_t l = 0;
00456 real_t sum = 0;
00457 for (size_t i=0; i<__basisX.dimension(); ++i) {
00458 const real_t vi = baseValuesX[i];
00459 for (size_t j=0;j<__basisY.dimension(); ++j) {
00460 const real_t vj = baseValuesY[j];
00461 const real_t vi_vj = vi*vj;
00462 for (size_t k=0; k<__basisZ.dimension(); ++k) {
00463 const real_t vk = baseDerivativeValuesZ[k];
00464 const real_t vi_vj_vk = vi_vj * vk;
00465 sum += vi_vj_vk*__values[l];
00466 l++;
00467 }
00468 }
00469 }
00470
00471 return sum * __transformZ.inverseDeterminant();
00472 }
00473
00479 SpectralFunction(ConstReferenceCounting<SpectralMesh> mesh)
00480
00481 : ScalarFunctionBase(ScalarFunctionBase::spectral),
00482 __values(mesh->numberOfVertices()),
00483 __mesh(mesh),
00484 __discretizationType(ScalarDiscretizationTypeBase::spectralLegendre),
00485 __intervalX(__mesh->shape().a()[0],__mesh->shape().b()[0]),
00486 __intervalY(__mesh->shape().a()[1],__mesh->shape().b()[1]),
00487 __intervalZ(__mesh->shape().a()[2],__mesh->shape().b()[2]),
00488 __transformX(__intervalX),
00489 __transformY(__intervalY),
00490 __transformZ(__intervalZ),
00491 __gaussLobattoX(GaussLobattoManager::instance().get(__mesh->degree(0)+1)),
00492 __gaussLobattoY(GaussLobattoManager::instance().get(__mesh->degree(1)+1)),
00493 __gaussLobattoZ(GaussLobattoManager::instance().get(__mesh->degree(2)+1)),
00494 __basisX(__mesh->degree(0)),
00495 __basisY(__mesh->degree(1)),
00496 __basisZ(__mesh->degree(2))
00497 {
00498 ;
00499 }
00500
00507 SpectralFunction(ConstReferenceCounting<SpectralMesh> mesh,
00508 const ScalarFunctionBase& f)
00509 : ScalarFunctionBase(ScalarFunctionBase::spectral),
00510 __values((mesh->degree(0)+1)*(mesh->degree(1)+1)*(mesh->degree(2)+1)),
00511 __mesh(mesh),
00512 __discretizationType(ScalarDiscretizationTypeBase::spectralLegendre),
00513 __intervalX(__mesh->shape().a()[0],__mesh->shape().b()[0]),
00514 __intervalY(__mesh->shape().a()[1],__mesh->shape().b()[1]),
00515 __intervalZ(__mesh->shape().a()[2],__mesh->shape().b()[2]),
00516 __transformX(__intervalX),
00517 __transformY(__intervalY),
00518 __transformZ(__intervalZ),
00519 __gaussLobattoX(GaussLobattoManager::instance().get(__mesh->degree(0)+1)),
00520 __gaussLobattoY(GaussLobattoManager::instance().get(__mesh->degree(1)+1)),
00521 __gaussLobattoZ(GaussLobattoManager::instance().get(__mesh->degree(2)+1)),
00522 __basisX(__mesh->degree(0)),
00523 __basisY(__mesh->degree(1)),
00524 __basisZ(__mesh->degree(2))
00525 {
00526 (*this) = f;
00527 }
00528
00535 explicit SpectralFunction(const SpectralFunction& f)
00536 : ScalarFunctionBase(f),
00537 __values(f.__values),
00538 __mesh(f.__mesh),
00539 __discretizationType(f.__discretizationType),
00540 __intervalX(f.__intervalX),
00541 __intervalY(f.__intervalY),
00542 __intervalZ(f.__intervalZ),
00543 __transformX(f.__transformX),
00544 __transformY(f.__transformY),
00545 __transformZ(f.__transformZ),
00546 __gaussLobattoX(f.__gaussLobattoX),
00547 __gaussLobattoY(f.__gaussLobattoY),
00548 __gaussLobattoZ(f.__gaussLobattoZ),
00549 __basisX(f.__basisX),
00550 __basisY(f.__basisY),
00551 __basisZ(f.__basisZ)
00552 {
00553 ;
00554 }
00555
00562 SpectralFunction(ConstReferenceCounting<SpectralMesh> mesh,
00563 const real_t& d)
00564 : ScalarFunctionBase(ScalarFunctionBase::spectral),
00565 __values((mesh->degree(0)+1)*(mesh->degree(1)+1)*(mesh->degree(2)+1)),
00566 __mesh(mesh),
00567 __discretizationType(ScalarDiscretizationTypeBase::spectralLegendre),
00568 __intervalX(__mesh->shape().a()[0],__mesh->shape().b()[0]),
00569 __intervalY(__mesh->shape().a()[1],__mesh->shape().b()[1]),
00570 __intervalZ(__mesh->shape().a()[2],__mesh->shape().b()[2]),
00571 __transformX(__intervalX),
00572 __transformY(__intervalY),
00573 __transformZ(__intervalZ),
00574 __gaussLobattoX(GaussLobattoManager::instance().get(__mesh->degree(0)+1)),
00575 __gaussLobattoY(GaussLobattoManager::instance().get(__mesh->degree(1)+1)),
00576 __gaussLobattoZ(GaussLobattoManager::instance().get(__mesh->degree(2)+1)),
00577 __basisX(__mesh->degree(0)),
00578 __basisY(__mesh->degree(1)),
00579 __basisZ(__mesh->degree(2))
00580 {
00581 (*this) = ScalarFunctionConstant(d);
00582 }
00583
00590 SpectralFunction(ConstReferenceCounting<SpectralMesh> mesh,
00591 const Vector<real_t>& values)
00592
00593 : ScalarFunctionBase(ScalarFunctionBase::spectral),
00594 __values((mesh->degree(0)+1)*(mesh->degree(1)+1)*(mesh->degree(2)+1)),
00595 __mesh(mesh),
00596 __discretizationType(ScalarDiscretizationTypeBase::spectralLegendre),
00597 __intervalX(__mesh->shape().a()[0],__mesh->shape().b()[0]),
00598 __intervalY(__mesh->shape().a()[1],__mesh->shape().b()[1]),
00599 __intervalZ(__mesh->shape().a()[2],__mesh->shape().b()[2]),
00600 __transformX(__intervalX),
00601 __transformY(__intervalY),
00602 __transformZ(__intervalZ),
00603 __gaussLobattoX(GaussLobattoManager::instance().get(__mesh->degree(0)+1)),
00604 __gaussLobattoY(GaussLobattoManager::instance().get(__mesh->degree(1)+1)),
00605 __gaussLobattoZ(GaussLobattoManager::instance().get(__mesh->degree(2)+1)),
00606 __basisX(__mesh->degree(0)),
00607 __basisY(__mesh->degree(1)),
00608 __basisZ(__mesh->degree(2))
00609 {
00610 ASSERT(__values.size() == values.size());
00611 __values = values;
00612 }
00613
00618 ~SpectralFunction()
00619 {
00620 ;
00621 }
00622 };
00623 #endif // SPECTRAL_FUNCTION_HPP
00624
