Added: Projection of general functions onto a specified basis.
Currently, only Greville point and continuous global L2 is supported. Skip (without error) the projection and norm calculation when the projection type is not supported.
This commit is contained in:
Knut Morten Okstad
parent
906abc41d4
commit
e45e120a18
|
|
@ -533,10 +533,17 @@ public:
|
|||
//! problems using internal integration point buffers (with history-dependent
|
||||
//! data, etc.) and that must be traversed in the same sequence each time.
|
||||
//! \note The implementation of this method is placed in GlbL2projector.C
|
||||
bool L2projection(Matrix& sField,
|
||||
const IntegrandBase& integrand,
|
||||
bool L2projection(Matrix& sField, IntegrandBase* integrand,
|
||||
const TimeDomain& time);
|
||||
|
||||
//! \brief Projects an explicit function using a continuous global L2-fit.
|
||||
//! \param[out] fVals Control point values of the function
|
||||
//! \param[in] function The function to project
|
||||
//! \param[in] t Current time
|
||||
//!
|
||||
//! \note The implementation of this method is placed in GlbL2projector.C
|
||||
bool L2projection(Matrix& fVals, FunctionBase* function, double t = 0.0);
|
||||
|
||||
|
||||
// Methods for result extraction
|
||||
// =============================
|
||||
|
|
|
|||
|
|
@ -17,6 +17,8 @@
|
|||
#include "FiniteElement.h"
|
||||
#include "TimeDomain.h"
|
||||
#include "ASMbase.h"
|
||||
#include "IntegrandBase.h"
|
||||
#include "Function.h"
|
||||
#include "Profiler.h"
|
||||
|
||||
|
||||
|
|
@ -39,31 +41,51 @@ public:
|
|||
|
||||
GlbL2& gl2Int; //!< The global L2-projection integrand
|
||||
LocalIntegral* elmData; //!< Element data associated with problem integrand
|
||||
IntVec mnpc; //!< Matrix of element nodal correspondance
|
||||
std::vector<size_t> elem_sizes; //!< Size of each basis on the element
|
||||
std::vector<size_t> basis_sizes; //!< Size of each basis on the patch
|
||||
|
||||
IntVec mnpc; //!< Matrix of element nodal correspondance
|
||||
uIntVec elem_sizes; //!< Size of each basis on the element
|
||||
uIntVec basis_sizes; //!< Size of each basis on the patch
|
||||
};
|
||||
|
||||
|
||||
GlbL2::GlbL2 (IntegrandBase& p, size_t n) : problem(p), A(SparseMatrix::SUPERLU)
|
||||
GlbL2::GlbL2 (IntegrandBase* p, size_t n) : A(SparseMatrix::SUPERLU)
|
||||
{
|
||||
problem = p;
|
||||
function = nullptr;
|
||||
|
||||
A.redim(n,n);
|
||||
B.redim(n*p.getNoFields(2));
|
||||
B.redim(n*p->getNoFields(2));
|
||||
}
|
||||
|
||||
|
||||
GlbL2::GlbL2 (FunctionBase* f, size_t n) : A(SparseMatrix::SUPERLU)
|
||||
{
|
||||
problem = nullptr;
|
||||
function = f;
|
||||
|
||||
A.redim(n,n);
|
||||
B.redim(n*f->dim());
|
||||
}
|
||||
|
||||
|
||||
int GlbL2::getIntegrandType () const
|
||||
{
|
||||
// Mask off the element interface flag
|
||||
return problem.getIntegrandType() & ~INTERFACE_TERMS;
|
||||
if (problem)
|
||||
// Mask off the element interface flag
|
||||
return problem->getIntegrandType() & ~INTERFACE_TERMS;
|
||||
else
|
||||
return STANDARD;
|
||||
}
|
||||
|
||||
|
||||
LocalIntegral* GlbL2::getLocalIntegral (size_t nen, size_t iEl,
|
||||
bool neumann) const
|
||||
{
|
||||
return new L2Mats(*const_cast<GlbL2*>(this),
|
||||
problem.getLocalIntegral(nen,iEl,neumann));
|
||||
if (problem)
|
||||
return new L2Mats(*const_cast<GlbL2*>(this),
|
||||
problem->getLocalIntegral(nen,iEl,neumann));
|
||||
else
|
||||
return new L2Mats(*const_cast<GlbL2*>(this));
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -74,13 +96,15 @@ bool GlbL2::initElement (const IntVec& MNPC, const FiniteElement& fe,
|
|||
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
|
||||
|
||||
gl2.mnpc = MNPC;
|
||||
return problem.initElement(MNPC,fe,Xc,nPt,*gl2.elmData);
|
||||
if (problem && gl2.elmData)
|
||||
return problem->initElement(MNPC,fe,Xc,nPt,*gl2.elmData);
|
||||
else
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool GlbL2::initElement (const IntVec& MNPC1,
|
||||
const std::vector<size_t>& elem_sizes,
|
||||
const std::vector<size_t>& basis_sizes,
|
||||
const uIntVec& elem_sizes, const uIntVec& basis_sizes,
|
||||
LocalIntegral& elmInt)
|
||||
{
|
||||
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
|
||||
|
|
@ -88,7 +112,10 @@ bool GlbL2::initElement (const IntVec& MNPC1,
|
|||
gl2.mnpc = MNPC1;
|
||||
gl2.elem_sizes = elem_sizes;
|
||||
gl2.basis_sizes = basis_sizes;
|
||||
return problem.initElement(MNPC1,elem_sizes,basis_sizes,*gl2.elmData);
|
||||
if (problem && gl2.elmData)
|
||||
return problem->initElement(MNPC1,elem_sizes,basis_sizes,*gl2.elmData);
|
||||
else
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -100,24 +127,16 @@ bool GlbL2::evalInt (LocalIntegral& elmInt,
|
|||
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
|
||||
|
||||
Vector solPt;
|
||||
if (!problem.evalSol(solPt,fe,X,gl2.mnpc))
|
||||
if (!problem.diverged(fe.iGP+1))
|
||||
return false;
|
||||
|
||||
size_t a, b, nnod = A.dim();
|
||||
for (a = 0; a < fe.N.size(); a++)
|
||||
if (problem)
|
||||
{
|
||||
int inod = gl2.mnpc[a]+1;
|
||||
for (b = 0; b < fe.N.size(); b++)
|
||||
{
|
||||
int jnod = gl2.mnpc[b]+1;
|
||||
A(inod,jnod) += fe.N[a]*fe.N[b]*fe.detJxW;
|
||||
}
|
||||
for (b = 0; b < solPt.size(); b++)
|
||||
B(inod+b*nnod) += fe.N[a]*solPt[b]*fe.detJxW;
|
||||
if (!problem->evalSol(solPt,fe,X,gl2.mnpc))
|
||||
if (!problem->diverged(fe.iGP+1))
|
||||
return false;
|
||||
}
|
||||
else if (function)
|
||||
solPt = function->getValue(X);
|
||||
|
||||
return true;
|
||||
return this->formL2Mats(gl2.mnpc,solPt,fe,X);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -129,21 +148,33 @@ bool GlbL2::evalIntMx (LocalIntegral& elmInt,
|
|||
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
|
||||
|
||||
Vector solPt;
|
||||
if (!problem.evalSol(solPt,fe,X,gl2.mnpc,gl2.elem_sizes,gl2.basis_sizes))
|
||||
if (!problem.diverged(fe.iGP+1))
|
||||
return false;
|
||||
|
||||
size_t a, b, nnod = A.dim();
|
||||
for (a = 0; a < fe.basis(1).size(); a++)
|
||||
if (problem)
|
||||
{
|
||||
int inod = gl2.mnpc[a]+1;
|
||||
for (b = 0; b < fe.basis(1).size(); b++)
|
||||
if (!problem->evalSol(solPt,fe,X,gl2.mnpc,gl2.elem_sizes,gl2.basis_sizes))
|
||||
if (!problem->diverged(fe.iGP+1))
|
||||
return false;
|
||||
}
|
||||
else if (function)
|
||||
solPt = function->getValue(X);
|
||||
|
||||
return this->formL2Mats(gl2.mnpc,solPt,fe,X);
|
||||
}
|
||||
|
||||
|
||||
bool GlbL2::formL2Mats (const IntVec& mnpc, const Vector& solPt,
|
||||
const FiniteElement& fe, const Vec3& X) const
|
||||
{
|
||||
size_t a, b, nnod = A.dim();
|
||||
for (a = 0; a < fe.N.size(); a++)
|
||||
{
|
||||
int inod = mnpc[a]+1;
|
||||
for (b = 0; b < fe.N.size(); b++)
|
||||
{
|
||||
int jnod = gl2.mnpc[b]+1;
|
||||
A(inod,jnod) += fe.basis(1)[a]*fe.basis(1)[b]*fe.detJxW;
|
||||
int jnod = mnpc[b]+1;
|
||||
A(inod,jnod) += fe.N[a]*fe.N[b]*fe.detJxW;
|
||||
}
|
||||
for (b = 0; b < solPt.size(); b++)
|
||||
B(inod+b*nnod) += fe.basis(1)[a]*solPt[b]*fe.detJxW;
|
||||
B(inod+b*nnod) += fe.N[a]*solPt[b]*fe.detJxW;
|
||||
}
|
||||
|
||||
return true;
|
||||
|
|
@ -173,7 +204,11 @@ bool GlbL2::solve (Matrix& sField)
|
|||
if (!A.solve(B)) return false;
|
||||
|
||||
// Store the nodal values of the projected field
|
||||
size_t j, ncomp = problem.getNoFields(2);
|
||||
size_t j, ncomp = 0;
|
||||
if (problem)
|
||||
ncomp = problem->getNoFields(2);
|
||||
else if (function)
|
||||
ncomp = function->dim();
|
||||
sField.resize(ncomp,nnod);
|
||||
for (i = 1; i <= nnod; i++)
|
||||
for (j = 1; j <= ncomp; j++)
|
||||
|
|
@ -187,12 +222,12 @@ bool GlbL2::solve (Matrix& sField)
|
|||
|
||||
|
||||
bool ASMbase::L2projection (Matrix& sField,
|
||||
const IntegrandBase& integrand,
|
||||
IntegrandBase* integrand,
|
||||
const TimeDomain& time)
|
||||
{
|
||||
PROFILE2("ASMbase::L2projection");
|
||||
|
||||
GlbL2 gl2(const_cast<IntegrandBase&>(integrand),this->getNoNodes(1));
|
||||
GlbL2 gl2(integrand,this->getNoNodes(1));
|
||||
GlobalIntegral dummy;
|
||||
|
||||
gl2.preAssemble(MNPC,this->getNoElms(true));
|
||||
|
|
@ -200,6 +235,19 @@ bool ASMbase::L2projection (Matrix& sField,
|
|||
}
|
||||
|
||||
|
||||
bool ASMbase::L2projection (Matrix& sField, FunctionBase* function, double t)
|
||||
{
|
||||
PROFILE2("ASMbase::L2projection");
|
||||
|
||||
GlbL2 gl2(function,this->getNoNodes(1));
|
||||
GlobalIntegral dummy;
|
||||
TimeDomain time; time.t = t;
|
||||
|
||||
gl2.preAssemble(MNPC,this->getNoElms(true));
|
||||
return this->integrate(gl2,dummy,time) && gl2.solve(sField);
|
||||
}
|
||||
|
||||
|
||||
bool ASMbase::globalL2projection (Matrix& sField,
|
||||
const IntegrandBase& integrand,
|
||||
bool continuous) const
|
||||
|
|
|
|||
|
|
@ -14,9 +14,14 @@
|
|||
#ifndef _GLB_L2_PROJECTOR_H
|
||||
#define _GLB_L2_PROJECTOR_H
|
||||
|
||||
#include "IntegrandBase.h"
|
||||
#include "Integrand.h"
|
||||
#include "SparseMatrix.h"
|
||||
|
||||
class IntegrandBase;
|
||||
class FunctionBase;
|
||||
|
||||
typedef std::vector<size_t> uIntVec; //!< General unsigned integer vector
|
||||
|
||||
|
||||
/*!
|
||||
\brief General integrand for L2-projection of secondary solutions.
|
||||
|
|
@ -28,7 +33,11 @@ public:
|
|||
//! \brief The constructor initializes the projection matrices.
|
||||
//! \param[in] p The main problem integrand
|
||||
//! \param[in] n Dimension of the L2-projection matrices (number of nodes)
|
||||
GlbL2(IntegrandBase& p, size_t n);
|
||||
GlbL2(IntegrandBase* p, size_t n);
|
||||
//! \brief Alternative constructor for projection of explicit functions.
|
||||
//! \param[in] f The function to to L2-projection on
|
||||
//! \param[in] n Dimension of the L2-projection matrices (number of nodes)
|
||||
GlbL2(FunctionBase* f, size_t n);
|
||||
//! \brief Empty destructor.
|
||||
virtual ~GlbL2() {}
|
||||
|
||||
|
|
@ -52,24 +61,22 @@ public:
|
|||
virtual bool initElement(const IntVec& MNPC, const FiniteElement& fe,
|
||||
const Vec3& X0, size_t nPt,
|
||||
LocalIntegral& elmInt);
|
||||
//! \brief Initializes current element for numerical integration (mixed integrands).
|
||||
//! \brief Initializes current element for numerical integration (mixed).
|
||||
//! \param[in] MNPC1 Matrix of nodal point correspondance for current element
|
||||
//! \param[in] elem_sizes Size of each basis on the element
|
||||
//! \param[in] basis_sizes Size of each basis on the patch
|
||||
//! \param elmInt Local integral for element
|
||||
virtual bool initElement(const IntVec& MNPC1,
|
||||
const std::vector<size_t>& elem_sizes,
|
||||
const std::vector<size_t>& basis_sizes,
|
||||
const uIntVec& elem_sizes,
|
||||
const uIntVec& basis_sizes,
|
||||
LocalIntegral& elmInt);
|
||||
|
||||
//! \brief Dummy implementation.
|
||||
virtual bool initElement(const IntVec&, LocalIntegral&) { return false; }
|
||||
|
||||
//! \brief Dummy implementation.
|
||||
virtual bool initElementBou(const IntVec&, LocalIntegral&) { return false; }
|
||||
//! \brief Dummy implementation.
|
||||
virtual bool initElementBou(const IntVec&, const std::vector<size_t>&,
|
||||
const std::vector<size_t>&,
|
||||
virtual bool initElementBou(const IntVec&, const uIntVec&, const uIntVec&,
|
||||
LocalIntegral&) { return false; }
|
||||
|
||||
using Integrand::evalInt;
|
||||
|
|
@ -97,8 +104,18 @@ public:
|
|||
//! \param[out] sField Nodal/control-point values of the projected results.
|
||||
bool solve(Matrix& sField);
|
||||
|
||||
protected:
|
||||
//! \brief Integrates the L2-projection matrices.
|
||||
//! \param[in] mnpc Matrix of nodal point correspondance
|
||||
//! \param[in] solPt Integration point values of the field to project
|
||||
//! \param[in] fe Finite element data of current integration point
|
||||
//! \param[in] X Cartesian coordinates of current integration point
|
||||
bool formL2Mats(const IntVec& mnpc, const Vector& solPt,
|
||||
const FiniteElement& fe, const Vec3& X) const;
|
||||
|
||||
private:
|
||||
IntegrandBase& problem; //!< The main problem integrand
|
||||
IntegrandBase* problem; //!< The main problem integrand
|
||||
FunctionBase* function; //!< Explicit function to L2-project
|
||||
mutable SparseMatrix A; //!< Left-hand-side matrix of the L2-projection
|
||||
mutable StdVector B; //!< Right-hand-side vectors of the L2-projection
|
||||
};
|
||||
|
|
|
|||
|
|
@ -66,8 +66,9 @@ SIMbase::~SIMbase ()
|
|||
IFEM::cout <<"\nEntering SIMbase destructor"<< std::endl;
|
||||
#endif
|
||||
|
||||
for (IntegrandMap::iterator it = myInts.begin(); it != myInts.end(); ++it)
|
||||
if (it->second != myProblem) delete it->second;
|
||||
for (auto& itg : myInts)
|
||||
if (itg.second != myProblem)
|
||||
delete itg.second;
|
||||
|
||||
if (myProblem) delete myProblem;
|
||||
if (mySol) delete mySol;
|
||||
|
|
@ -75,8 +76,8 @@ SIMbase::~SIMbase ()
|
|||
if (mySam) delete mySam;
|
||||
if (mySolParams) delete mySolParams;
|
||||
|
||||
for (PatchVec::iterator i1 = myModel.begin(); i1 != myModel.end(); i1++)
|
||||
delete *i1;
|
||||
for (ASMbase* patch : myModel)
|
||||
delete patch;
|
||||
|
||||
myModel.clear();
|
||||
this->SIMbase::clearProperties();
|
||||
|
|
@ -89,14 +90,15 @@ SIMbase::~SIMbase ()
|
|||
|
||||
void SIMbase::clearProperties ()
|
||||
{
|
||||
for (PatchVec::iterator i1 = myModel.begin(); i1 != myModel.end(); i1++)
|
||||
(*i1)->clear(true); // retain the geometry only
|
||||
for (SclFuncMap::iterator i2 = myScalars.begin(); i2 != myScalars.end(); i2++)
|
||||
delete i2->second;
|
||||
for (VecFuncMap::iterator i3 = myVectors.begin(); i3 != myVectors.end(); i3++)
|
||||
delete i3->second;
|
||||
for (TracFuncMap::iterator i4 = myTracs.begin(); i4 != myTracs.end(); i4++)
|
||||
delete i4->second;
|
||||
for (ASMbase* patch : myModel)
|
||||
patch->clear(true); // retain the geometry only
|
||||
|
||||
for (auto& i2 : myScalars)
|
||||
delete i2.second;
|
||||
for (auto& i3 : myVectors)
|
||||
delete i3.second;
|
||||
for (auto& i4 : myTracs)
|
||||
delete i4.second;
|
||||
|
||||
myPatches.clear();
|
||||
myGlb2Loc.clear();
|
||||
|
|
@ -164,13 +166,12 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
|
|||
if (!this->createFEMmodel('Y')) return false;
|
||||
|
||||
PatchVec::const_iterator mit;
|
||||
IntVec::const_iterator it;
|
||||
ASMbase* pch;
|
||||
size_t patch;
|
||||
|
||||
// Erase all patches that should be ignored in the analysis
|
||||
for (it = ignored.begin(); it != ignored.end(); ++it)
|
||||
if ((pch = this->getPatch(*it)))
|
||||
for (int idx : ignored)
|
||||
if ((pch = this->getPatch(idx)))
|
||||
pch->clear();
|
||||
|
||||
// If material properties are specified for at least one patch, assign the
|
||||
|
|
@ -1257,7 +1258,9 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
|||
norm->initIntegration(nIntGP,nBouGP);
|
||||
|
||||
// Number of recovered solution components
|
||||
size_t nCmp = ssol.empty() ? 0 : ssol.front().size() / mySam->getNoNodes();
|
||||
size_t i, nCmp = 0;
|
||||
for (i = 0; i < ssol.size() && nCmp == 0; i++)
|
||||
nCmp = ssol[i].size() / mySam->getNoNodes();
|
||||
|
||||
#ifdef USE_OPENMP
|
||||
// When assembling in parallel, we must always do the norm summation
|
||||
|
|
@ -1270,11 +1273,13 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
|||
// Initialize norm integral classes
|
||||
gNorm.resize(norm->getNoFields(0));
|
||||
size_t nNorms = 0;
|
||||
for (size_t j = 0; j < gNorm.size(); j++) {
|
||||
size_t nNrm = norm->getNoFields(1+j);
|
||||
gNorm[j].resize(nNrm,true);
|
||||
nNorms += nNrm;
|
||||
}
|
||||
for (i = 0; i < gNorm.size(); i++)
|
||||
if (i == 0 || i > ssol.size() || !ssol[i-1].empty())
|
||||
{
|
||||
size_t nNrm = norm->getNoFields(1+i);
|
||||
gNorm[i].resize(nNrm,true);
|
||||
nNorms += nNrm;
|
||||
}
|
||||
|
||||
GlbNorm globalNorm(gNorm,norm->getFinalOperation());
|
||||
LintegralVec elementNorms;
|
||||
|
|
@ -1282,7 +1287,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
|||
{
|
||||
eNorm->resize(nNorms,mySam->getNoElms(),true);
|
||||
elementNorms.reserve(eNorm->cols());
|
||||
for (size_t i = 0; i < eNorm->cols(); i++)
|
||||
for (i = 0; i < eNorm->cols(); i++)
|
||||
elementNorms.push_back(new ElmNorm(eNorm->ptr(i),nNorms));
|
||||
norm->setLocalIntegrals(&elementNorms);
|
||||
}
|
||||
|
|
@ -1302,8 +1307,10 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
|||
lp = p->patch;
|
||||
ok = this->extractPatchSolution(psol,lp-1);
|
||||
for (k = 0; k < ssol.size(); k++)
|
||||
if (!ssol[k].empty())
|
||||
pch->extractNodeVec(ssol[k],norm->getProjection(k),nCmp);
|
||||
if (ssol[k].empty())
|
||||
norm->getProjection(k).clear();
|
||||
else
|
||||
pch->extractNodeVec(ssol[k],norm->getProjection(k),nCmp);
|
||||
ok &= pch->integrate(*norm,globalNorm,time);
|
||||
}
|
||||
else
|
||||
|
|
@ -1312,12 +1319,14 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
|||
if (lp == 0)
|
||||
// All patches are referring to the same material, and we assume it has
|
||||
// been initialized during input processing (thus no initMaterial call here)
|
||||
for (size_t i = 0; i < myModel.size() && ok; i++)
|
||||
for (i = 0; i < myModel.size() && ok; i++)
|
||||
{
|
||||
ok = this->extractPatchSolution(psol,i);
|
||||
for (k = 0; k < ssol.size(); k++)
|
||||
if (!ssol[k].empty())
|
||||
myModel[i]->extractNodeVec(ssol[k],norm->getProjection(k),nCmp);
|
||||
if (ssol[k].empty())
|
||||
norm->getProjection(k).clear();
|
||||
else
|
||||
myModel[i]->extractNodeVec(ssol[k],norm->getProjection(k),nCmp);
|
||||
ok &= myModel[i]->integrate(*norm,globalNorm,time);
|
||||
lp = i+1;
|
||||
}
|
||||
|
|
@ -1560,7 +1569,7 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
|
|||
case SIMoptions::CGL2_INT:
|
||||
if (msgLevel > 1 && i == 0)
|
||||
IFEM::cout <<"\tContinuous global L2-projection"<< std::endl;
|
||||
ok = myModel[i]->L2projection(values,*myProblem,time);
|
||||
ok = myModel[i]->L2projection(values,myProblem,time);
|
||||
break;
|
||||
|
||||
case SIMoptions::SCR:
|
||||
|
|
@ -1628,27 +1637,52 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
|
|||
}
|
||||
|
||||
|
||||
bool SIMbase::project (Vector& values, const RealFunc* f,
|
||||
int basis, int iField, int nFields, double time) const
|
||||
bool SIMbase::project (Vector& values, const FunctionBase* f,
|
||||
int basis, int iField, int nFields,
|
||||
SIMoptions::ProjectionMethod pMethod, double time) const
|
||||
{
|
||||
bool ok = true;
|
||||
for (size_t j = 0; j < myModel.size() && ok; j++)
|
||||
{
|
||||
if (myModel[j]->empty()) continue; // skip empty patches
|
||||
|
||||
Vector loc_scalar;
|
||||
ok = myModel[j]->evaluate(f,loc_scalar,basis,time);
|
||||
Vector loc_values;
|
||||
Matrix f_values;
|
||||
switch (pMethod) {
|
||||
case SIMoptions::GLOBAL:
|
||||
// Greville point projection
|
||||
ok = myModel[j]->evaluate(f,loc_values,basis,time);
|
||||
break;
|
||||
|
||||
if (nFields <= 1)
|
||||
ok &= myModel[j]->injectNodeVec(loc_scalar,values,1,basis);
|
||||
case SIMoptions::CGL2:
|
||||
case SIMoptions::CGL2_INT:
|
||||
// Continuous global L2-projection
|
||||
ok = myModel[j]->L2projection(f_values,const_cast<FunctionBase*>(f),time);
|
||||
loc_values = f_values;
|
||||
break;
|
||||
|
||||
default:
|
||||
std::cerr <<" *** SIMbase::project: Projection method "<< pMethod
|
||||
<<" not implemented for functions."<< std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
if (nFields <= (int)f->dim())
|
||||
ok &= myModel[j]->injectNodeVec(loc_values,values,f->dim(),basis);
|
||||
else if (f->dim() > 1)
|
||||
{
|
||||
std::cerr <<" *** SIMbase::project: Cannot interleave non-scalar function"
|
||||
<<" into a multi-dimensional field."<< std::endl;
|
||||
return false;
|
||||
}
|
||||
else if (iField < nFields)
|
||||
{
|
||||
// Interleave
|
||||
size_t i, k = iField;
|
||||
Vector loc_vector(loc_scalar.size()*nFields);
|
||||
Vector loc_vector(loc_values.size()*nFields);
|
||||
myModel[j]->extractNodeVec(values,loc_vector,0,basis);
|
||||
for (i = 0; i < loc_scalar.size(); i++, k += nFields)
|
||||
loc_vector[k] = loc_scalar[i];
|
||||
for (i = 0; i < loc_values.size(); i++, k += nFields)
|
||||
loc_vector[k] = loc_values[i];
|
||||
ok &= myModel[j]->injectNodeVec(loc_vector,values,0,basis);
|
||||
}
|
||||
else
|
||||
|
|
|
|||
|
|
@ -28,6 +28,7 @@ class SAM;
|
|||
class AlgEqSystem;
|
||||
class LinSolParams;
|
||||
class SystemVector;
|
||||
class FunctionBase;
|
||||
class RealFunc;
|
||||
class VecFunc;
|
||||
class TractionFunc;
|
||||
|
|
@ -419,15 +420,17 @@ public:
|
|||
bool project(Vector& ssol, const Vector& psol,
|
||||
SIMoptions::ProjectionMethod pMethod = SIMoptions::GLOBAL) const;
|
||||
|
||||
//! \brief Projects a scalar function onto the specified basis.
|
||||
//! \brief Projects a function onto the specified basis.
|
||||
//! \param[out] values Resulting control point values
|
||||
//! \param[in] f The function to evaluate
|
||||
//! \param[in] basis Which basis to consider
|
||||
//! \param[in] iField Field component offset in values vector
|
||||
//! \param[in] nFields Number of field components in values vector
|
||||
//! \param[in] pMethod Projection method to use
|
||||
//! \param[in] time Current time
|
||||
bool project(Vector& values, const RealFunc* f,
|
||||
bool project(Vector& values, const FunctionBase* f,
|
||||
int basis = 1, int iField = 0, int nFields = 1,
|
||||
SIMoptions::ProjectionMethod pMethod = SIMoptions::GLOBAL,
|
||||
double time = 0.0) const;
|
||||
|
||||
//! \brief Evaluates the secondary solution field for specified patch.
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user