Added calculation of norms for projected secondary solution, etc.
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1120 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
Knut Morten Okstad
parent
d798dfcd6c
commit
ebd0512817
@ -233,9 +233,9 @@ int main (int argc, char** argv)
|
||||
|
||||
model->setQuadratureRule(nGauss);
|
||||
|
||||
Matrix eNorm;
|
||||
Matrix eNorm, ssol;
|
||||
Vector gNorm, load;
|
||||
Vectors displ(1);
|
||||
Vectors displ(1), projs;
|
||||
std::vector<Mode> modes;
|
||||
std::vector<Mode>::const_iterator it;
|
||||
|
||||
@ -255,20 +255,40 @@ int main (int argc, char** argv)
|
||||
if (!model->solveSystem(displ.front(),1))
|
||||
return 3;
|
||||
|
||||
// Evaluate solution norms
|
||||
model->setMode(SIM::RECOVERY);
|
||||
if (!model->solutionNorms(displ,eNorm,gNorm))
|
||||
if (SIMbase::discretization == SIMbase::Spline)
|
||||
{
|
||||
// Project the FE stresses onto the splines basis
|
||||
if (!model->project(ssol,displ.front(),SIMbase::GLOBAL))
|
||||
return 4;
|
||||
else
|
||||
projs.push_back(ssol);
|
||||
}
|
||||
|
||||
// Integrate solution norms and errors
|
||||
if (!model->solutionNorms(displ,projs,eNorm,gNorm))
|
||||
return 4;
|
||||
|
||||
if (linalg.myPid == 0)
|
||||
{
|
||||
std::cout <<"Energy norm |u^h| = a(u^h,u^h)^0.5 : "<< gNorm(1);
|
||||
std::cout <<"\nExternal energy ((f,u^h)+(t,u^h)^0.5 : "<< gNorm(2);
|
||||
if (gNorm.size() > 2)
|
||||
std::cout <<"\nExact norm |u| = a(u,u)^0.5 : "<< gNorm(3);
|
||||
if (gNorm.size() > 3)
|
||||
std::cout <<"\nExact error a(e,e)^0.5, e=u-u^h : "<< gNorm(4)
|
||||
std::cout <<"\nExternal energy ((f,u^h)+(t,u^h)^0.5 : "<< gNorm(2);
|
||||
if (model->haveAnaSol() && gNorm.size() >= 4)
|
||||
std::cout <<"\nExact norm |u| = a(u,u)^0.5 : "<< gNorm(3)
|
||||
<<"\nExact error a(e,e)^0.5, e=u-u^h : "<< gNorm(4)
|
||||
<<"\nExact relative error (%) : "<< gNorm(4)/gNorm(3)*100.0;
|
||||
size_t j = model->haveAnaSol() ? 5 : 3;
|
||||
for (size_t i = 0; i < projs.size() && j < gNorm.size(); i++)
|
||||
{
|
||||
std::cout <<"\nEnergy norm |u^r| = a(u^r,u^r)^0.5 : "<< gNorm(j++);
|
||||
std::cout <<"\nError norm a(e,e)^0.5, e=u^r-u^h : "<< gNorm(j++);
|
||||
std::cout <<"\n- relative error (% of |u^r|) : "
|
||||
<< gNorm(j-1)/gNorm(j-2)*100.0;
|
||||
if (model->haveAnaSol() && j++ <= gNorm.size())
|
||||
std::cout <<"\nExact error a(e,e)^0.5, e=u-u^r : "<< gNorm(j-1)
|
||||
<<"\n- relative error (% of |u|) : "
|
||||
<< gNorm(j-1)/gNorm(3)*100.0;
|
||||
}
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
@ -349,6 +369,10 @@ int main (int argc, char** argv)
|
||||
if (!model->writeGlvS(displ.front(),n,iStep,nBlock))
|
||||
return 10;
|
||||
|
||||
// Write projected solution fields to VTF-file
|
||||
if (!model->writeGlvP(projs.front(),n,iStep,nBlock))
|
||||
return 10;
|
||||
|
||||
// Write eigenmodes
|
||||
for (it = modes.begin(); it != modes.end(); it++)
|
||||
if (!model->writeGlvM(*it, iop==3 || iop==4 || iop==6, n, nBlock))
|
||||
@ -359,7 +383,7 @@ int main (int argc, char** argv)
|
||||
if (!model->writeGlvN(eNorm,iStep,nBlock))
|
||||
return 12;
|
||||
|
||||
model->writeGlvStep(1,0,1);
|
||||
model->writeGlvStep(1);
|
||||
model->closeGlv();
|
||||
}
|
||||
|
||||
|
||||
@ -28,7 +28,7 @@ class ElmNorm : public LocalIntegral
|
||||
{
|
||||
public:
|
||||
//! \brief The constructor assigns the internal pointer.
|
||||
ElmNorm(double* p) : ptr(p) {}
|
||||
ElmNorm(double* p, size_t n) : ptr(p), nnv(n) {}
|
||||
//! \brief Empty destructor.
|
||||
virtual ~ElmNorm() {}
|
||||
|
||||
@ -37,8 +37,12 @@ public:
|
||||
//! \brief Indexing operator for referencing.
|
||||
const double& operator[](size_t i) const { return ptr[i]; }
|
||||
|
||||
//! \brief Returns the number of norm values.
|
||||
size_t size() const { return nnv; }
|
||||
|
||||
private:
|
||||
double* ptr; //!< Pointer to the actual norm values
|
||||
size_t nnv; //!< Number of norm values
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -35,9 +35,10 @@ bool GlbNorm::assemble (const LocalIntegral* elmObj, int elmId)
|
||||
if (!ptr) return false;
|
||||
|
||||
ElmNorm& elVals = *const_cast<ElmNorm*>(ptr);
|
||||
for (size_t i = 0; i < myVals.size(); i++)
|
||||
for (size_t i = 0; i < elVals.size(); i++)
|
||||
{
|
||||
myVals[i] += elVals[i];
|
||||
if (i < myVals.size())
|
||||
myVals[i] += elVals[i];
|
||||
this->applyFinalOp(elVals[i]);
|
||||
}
|
||||
|
||||
|
||||
@ -185,8 +185,8 @@ ElmNorm& NormBase::getElmNormBuffer (LocalIntegral*& elmInt, const size_t nn)
|
||||
|
||||
static double* data = 0;
|
||||
if (!data) data = new double[nn];
|
||||
static ElmNorm buf(data);
|
||||
memset(data,0,nn*sizeof(double));
|
||||
static ElmNorm buf(data,nn);
|
||||
elmInt = &buf;
|
||||
return buf;
|
||||
}
|
||||
|
||||
@ -221,6 +221,9 @@ public:
|
||||
//! \brief Returns the number of field components.
|
||||
virtual size_t getNoFields() const { return 0; }
|
||||
|
||||
//! \brief Accesses a projected secondary solution vector of current patch.
|
||||
virtual Vector& getProjection(size_t) { static Vector dummy; return dummy; }
|
||||
|
||||
protected:
|
||||
//! \brief Returns the element norm object to use in the integration.
|
||||
//! \param elmInt The local integral object to receive norm contributions
|
||||
|
||||
@ -550,9 +550,49 @@ NormBase* Elasticity::getNormIntegrand (AnaSol* asol) const
|
||||
}
|
||||
|
||||
|
||||
ElasticityNorm::ElasticityNorm (Elasticity& p, STensorFunc* a)
|
||||
: NormBase(p), anasol(a)
|
||||
{
|
||||
nrcmp = myProblem.getNoFields(2);
|
||||
}
|
||||
|
||||
|
||||
size_t ElasticityNorm::getNoFields () const
|
||||
{
|
||||
return anasol ? 4 : 2;
|
||||
size_t nf = anasol ? 4 : 2;
|
||||
for (size_t i = 0; i < prjsol.size(); i++)
|
||||
if (!prjsol.empty())
|
||||
nf += anasol ? 3 : 2;
|
||||
|
||||
return nf;
|
||||
}
|
||||
|
||||
|
||||
Vector& ElasticityNorm::getProjection (size_t i)
|
||||
{
|
||||
if (prjsol.size() < i)
|
||||
{
|
||||
prjsol.resize(i);
|
||||
mySols.resize(i);
|
||||
}
|
||||
|
||||
return prjsol[i-1];
|
||||
}
|
||||
|
||||
|
||||
bool ElasticityNorm::initElement (const std::vector<int>& MNPC)
|
||||
{
|
||||
// Extract projected solution vectors for this element
|
||||
int ierr = 0;
|
||||
for (size_t i = 0; i < mySols.size() && ierr == 0; i++)
|
||||
if (!prjsol[i].empty())
|
||||
ierr = utl::gather(MNPC,nrcmp,prjsol[i],mySols[i]);
|
||||
|
||||
if (ierr == 0) return myProblem.initElement(MNPC);
|
||||
|
||||
std::cerr <<" *** ElasticityNorm::initElement: Detected "
|
||||
<< ierr <<" node numbers out of range."<< std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
@ -560,18 +600,19 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
|
||||
const Vec3& X) const
|
||||
{
|
||||
Elasticity& problem = static_cast<Elasticity&>(myProblem);
|
||||
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
|
||||
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt,this->getNoFields());
|
||||
|
||||
// Evaluate the inverse constitutive matrix at this point
|
||||
Matrix Cinv;
|
||||
if (!problem.formCinverse(Cinv,X)) return false;
|
||||
|
||||
// Evaluate the finite element stress field
|
||||
Vector sigmah, sigma;
|
||||
Vector sigmah, sigma, error;
|
||||
if (!problem.evalSol(sigmah,fe.N,fe.dNdX,X))
|
||||
return false;
|
||||
else if (sigmah.size() == 4 && Cinv.rows() == 3)
|
||||
sigmah.erase(sigmah.begin()+2); // Remove the sigma_zz if plane strain
|
||||
|
||||
bool planeStrain = sigmah.size() == 4 && Cinv.rows() == 3;
|
||||
if (planeStrain) sigmah.erase(sigmah.begin()+2); // Remove the sigma_zz
|
||||
|
||||
double detJW = fe.detJxW;
|
||||
if (problem.isAxiSymmetric())
|
||||
@ -590,6 +631,7 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
|
||||
pnorm[1] += f*u*detJW;
|
||||
}
|
||||
|
||||
size_t ip = 2;
|
||||
if (anasol)
|
||||
{
|
||||
// Evaluate the analytical stress field
|
||||
@ -598,12 +640,36 @@ bool ElasticityNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
|
||||
sigma.erase(sigma.begin()+2); // Remove the sigma_zz if plane strain
|
||||
|
||||
// Integrate the energy norm a(u,u)
|
||||
pnorm[2] += sigma.dot(Cinv*sigma)*detJW;
|
||||
pnorm[ip++] += sigma.dot(Cinv*sigma)*detJW;
|
||||
// Integrate the error in energy norm a(u-u^h,u-u^h)
|
||||
sigma -= sigmah;
|
||||
pnorm[3] += sigma.dot(Cinv*sigma)*detJW;
|
||||
error = sigma - sigmah;
|
||||
pnorm[ip++] += error.dot(Cinv*error)*detJW;
|
||||
}
|
||||
|
||||
size_t i, j, k;
|
||||
for (i = 0; i < mySols.size(); i++)
|
||||
if (!mySols[i].empty())
|
||||
{
|
||||
// Evaluate projected stress field
|
||||
Vector sigmar(sigmah.size());
|
||||
for (j = k = 0; j < nrcmp && k < sigmar.size(); j++)
|
||||
if (!planeStrain || j != 2)
|
||||
sigmar[k++] = mySols[i].dot(fe.N,j,nrcmp);
|
||||
|
||||
// Integrate the energy norm a(u^r,u^r)
|
||||
pnorm[ip++] += sigmar.dot(Cinv*sigmar)*detJW;
|
||||
// Integrate the error in energy norm a(u^r-u^h,u^r-u^h)
|
||||
error = sigmar - sigmah;
|
||||
pnorm[ip++] += error.dot(Cinv*error)*detJW;
|
||||
|
||||
if (anasol)
|
||||
{
|
||||
// Integrate the error in the projected solution a(u-u^r,u-u^r)
|
||||
error = sigma - sigmar;
|
||||
pnorm[ip++] += error.dot(Cinv*error)*detJW;
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
@ -225,7 +225,7 @@ public:
|
||||
//! \brief The only constructor initializes its data members.
|
||||
//! \param[in] p The linear elasticity problem to evaluate norms for
|
||||
//! \param[in] a The analytical stress field (optional)
|
||||
ElasticityNorm(Elasticity& p, STensorFunc* a = 0) : NormBase(p), anasol(a) {}
|
||||
ElasticityNorm(Elasticity& p, STensorFunc* a = 0);
|
||||
//! \brief Empty destructor.
|
||||
virtual ~ElasticityNorm() {}
|
||||
|
||||
@ -234,6 +234,14 @@ public:
|
||||
//! \brief Returns a 1-based index of the external energy norm.
|
||||
virtual size_t indExt() const { return 2; }
|
||||
|
||||
//! \brief Initializes current element for numerical integration.
|
||||
//! \param[in] MNPC Matrix of nodal point correspondance for current element
|
||||
virtual bool initElement(const std::vector<int>& MNPC);
|
||||
//! \brief Initializes current element for numerical integration.
|
||||
//! \param[in] MNPC Matrix of nodal point correspondance for current element
|
||||
virtual bool initElement(const std::vector<int>& MNPC, const Vec3&, size_t)
|
||||
{ return this->initElement(MNPC); }
|
||||
|
||||
//! \brief Evaluates the integrand at an interior point.
|
||||
//! \param elmInt The local integral object to receive the contributions
|
||||
//! \param[in] fe Finite element data of current integration point
|
||||
@ -252,8 +260,16 @@ public:
|
||||
//! \brief Returns the number of norm quantities.
|
||||
virtual size_t getNoFields() const;
|
||||
|
||||
//! \brief Accesses a projected secondary solution vector of current patch.
|
||||
virtual Vector& getProjection(size_t i);
|
||||
|
||||
private:
|
||||
STensorFunc* anasol; //!< Analytical stress field
|
||||
|
||||
Vectors prjsol; //!< Projected secondary solution vectors for current patch
|
||||
Vectors mySols; //!< Local element vectors
|
||||
|
||||
size_t nrcmp; //!< Number of projected solution components
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -697,7 +697,8 @@ double SIMbase::solutionNorms (const Vector& x, double* inf,
|
||||
}
|
||||
|
||||
|
||||
bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
|
||||
bool SIMbase::solutionNorms (const TimeDomain& time,
|
||||
const Vectors& psol, const Vectors& ssol,
|
||||
Vector& gNorm, Matrix* eNorm)
|
||||
{
|
||||
NormBase* norm = myProblem->getNormIntegrand(mySol);
|
||||
@ -714,6 +715,12 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
|
||||
myProblem->initIntegration(time);
|
||||
const Vector& primsol = psol.front();
|
||||
|
||||
size_t i, j, k;
|
||||
for (i = 0; i < ssol.size(); i++)
|
||||
if (!ssol[i].empty())
|
||||
norm->getProjection(i+1);
|
||||
|
||||
size_t nCmp = ssol.empty() ? 0 : ssol.front().size() / mySam->getNoNodes();
|
||||
size_t nNorms = norm->getNoFields();
|
||||
gNorm.resize(nNorms,true);
|
||||
|
||||
@ -724,21 +731,24 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
|
||||
{
|
||||
eNorm->resize(nNorms,mySam->getNoElms(),true);
|
||||
elementNorms.reserve(eNorm->cols());
|
||||
for (size_t i = 0; i < eNorm->cols(); i++)
|
||||
elementNorms.push_back(new ElmNorm(eNorm->ptr(i)));
|
||||
for (i = 0; i < eNorm->cols(); i++)
|
||||
elementNorms.push_back(new ElmNorm(eNorm->ptr(i),nNorms));
|
||||
}
|
||||
|
||||
// Loop over the different material regions, integrating solution norm terms
|
||||
// for the patch domain associated with each material
|
||||
bool ok = true;
|
||||
size_t i, j = 0, lp = 0;
|
||||
for (i = 0; i < myProps.size() && ok; i++)
|
||||
size_t lp = 0;
|
||||
for (i = j = 0; i < myProps.size() && ok; i++)
|
||||
if (myProps[i].pcode == Property::MATERIAL)
|
||||
if ((j = myProps[i].patch) < 1 || j > myModel.size())
|
||||
ok = false;
|
||||
else if (this->initMaterial(myProps[i].pindx))
|
||||
{
|
||||
myModel[j-1]->extractNodeVec(primsol,myProblem->getSolution());
|
||||
for (k = 0; k < ssol.size(); k++)
|
||||
if (!ssol[k].empty())
|
||||
myModel[j-1]->extractNodeVec(ssol[k],norm->getProjection(k+1),nCmp);
|
||||
ok = myModel[j-1]->integrate(*norm,globalNorm,time,elementNorms);
|
||||
lp = j;
|
||||
}
|
||||
@ -751,8 +761,11 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
|
||||
for (i = 0; i < myModel.size() && ok; i++)
|
||||
{
|
||||
myModel[i]->extractNodeVec(primsol,myProblem->getSolution());
|
||||
for (k = 0; k < ssol.size(); k++)
|
||||
if (!ssol[k].empty())
|
||||
myModel[i]->extractNodeVec(ssol[k],norm->getProjection(k+1),nCmp);
|
||||
ok = myModel[i]->integrate(*norm,globalNorm,time,elementNorms);
|
||||
lp = j;
|
||||
lp = i+1;
|
||||
}
|
||||
|
||||
// Integrate norm contributions due to Neumann boundary conditions, if any.
|
||||
@ -900,13 +913,15 @@ bool SIMbase::writeGlv (const char* inpFile, const int* nViz, int format)
|
||||
{
|
||||
if (myVtf) return false;
|
||||
|
||||
#if HAS_VTFAPI == 2
|
||||
const char* ext = ".vtfx";
|
||||
#else
|
||||
const char* ext = ".vtf";
|
||||
#endif
|
||||
|
||||
// Open a new VTF-file
|
||||
char* vtfName = new char[strlen(inpFile)+10];
|
||||
strtok(strcpy(vtfName,inpFile),".");
|
||||
const char* ext = ".vtf";
|
||||
#if HAS_VTFAPI == 2
|
||||
ext = ".vtfx";
|
||||
#endif
|
||||
if (nProc > 1)
|
||||
sprintf(vtfName+strlen(vtfName),"_p%04d%s",myPid,ext);
|
||||
else
|
||||
@ -1220,6 +1235,52 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
|
||||
}
|
||||
|
||||
|
||||
bool SIMbase::writeGlvP (const Vector& ssol, const int* nViz,
|
||||
int iStep, int& nBlock, double time,
|
||||
int idBlock, const char* prefix)
|
||||
{
|
||||
if (ssol.empty())
|
||||
return true;
|
||||
else if (!myVtf)
|
||||
return false;
|
||||
|
||||
Matrix field;
|
||||
Vector lovec;
|
||||
const size_t nf = myProblem->getNoFields(2);
|
||||
std::vector<int> sID[nf];
|
||||
|
||||
size_t i, j;
|
||||
int geomID = 0;
|
||||
for (i = 0; i < myModel.size(); i++)
|
||||
{
|
||||
if (myModel[i]->empty()) continue; // skip empty patches
|
||||
|
||||
if (msgLevel > 1)
|
||||
std::cout <<"Writing projected solution for patch "<< i+1 << std::endl;
|
||||
|
||||
// Evaluate the solution variables at the visualization points
|
||||
myModel[i]->extractNodeVec(ssol,lovec,nf);
|
||||
if (!myModel[i]->evalSolution(field,lovec,nViz))
|
||||
return false;
|
||||
|
||||
// Write out to VTF-file as scalar fields
|
||||
geomID++;
|
||||
for (j = 0; j < field.rows(); j++)
|
||||
if (!myVtf->writeNres(field.getRow(1+j),++nBlock,geomID))
|
||||
return false;
|
||||
else
|
||||
sID[j].push_back(nBlock);
|
||||
}
|
||||
|
||||
// Write result block identifications
|
||||
for (j = 0; j < nf && !sID[j].empty(); j++)
|
||||
if (!myVtf->writeSblk(sID[j],myProblem->getField2Name(j,prefix),
|
||||
++idBlock,iStep)) return false;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool SIMbase::writeGlvStep (int iStep, double value, int itype)
|
||||
{
|
||||
if (itype == 0)
|
||||
@ -1280,7 +1341,7 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
|
||||
|
||||
Matrix field;
|
||||
int geomID = 0;
|
||||
std::vector<int> sID[3];
|
||||
std::vector<int> sID[10];
|
||||
|
||||
size_t i, j, k;
|
||||
for (i = 0; i < myModel.size(); i++)
|
||||
@ -1293,25 +1354,39 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock)
|
||||
geomID++;
|
||||
myModel[i]->extractElmRes(norms,field);
|
||||
|
||||
for (j = k = 0; j < field.rows() && j < 4; j++)
|
||||
if (field.rows()%2 || j != 1) // Skip the external norms (always zero)
|
||||
for (j = k = 0; j < field.rows() && j < 10; j++)
|
||||
if (j != 1) // Skip the external norms (always zero)
|
||||
if (!myVtf->writeEres(field.getRow(1+j),++nBlock,geomID))
|
||||
return false;
|
||||
else
|
||||
sID[k++].push_back(nBlock);
|
||||
}
|
||||
|
||||
const char* label[3] = {
|
||||
const char* label[9] = {
|
||||
"a(u^h,u^h)^0.5",
|
||||
"a(u,u)^0.5",
|
||||
"a(e,e)^0.5, e=u-u^h"
|
||||
"a(e,e)^0.5, e=u-u^h",
|
||||
"a(u^r,u^r)^0.5",
|
||||
"a(e,e)^0.5, e=u^r-u^h",
|
||||
"a(e,e)^0.5, e=u-u^r",
|
||||
"a(u^rr,u^rr)^0.5",
|
||||
"a(e,e)^0.5, e=u^rr-u^h",
|
||||
"a(e,e)^0.5, e=u-u^rr"
|
||||
};
|
||||
|
||||
int idBlock = 100;
|
||||
for (j = 0; j < 3; j++)
|
||||
if (!sID[j].empty())
|
||||
if (!myVtf->writeSblk(sID[j],label[j],++idBlock,iStep))
|
||||
return false;
|
||||
int idBlock = 200;
|
||||
for (j = k = 0; !sID[j].empty(); j++, k++)
|
||||
{
|
||||
if (!mySol)
|
||||
{
|
||||
if (k == 1)
|
||||
k = 3;
|
||||
else if (k%3 == 2)
|
||||
k ++;
|
||||
}
|
||||
if (!myVtf->writeSblk(sID[j],label[k],++idBlock,iStep))
|
||||
return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
@ -1375,9 +1450,9 @@ bool SIMbase::dumpGeometry (std::ostream& os) const
|
||||
|
||||
bool SIMbase::dumpBasis (std::ostream& os, int basis, size_t patch) const
|
||||
{
|
||||
size_t start = patch?patch-1:0;
|
||||
size_t end = patch?start+1:myModel.size();
|
||||
for (size_t i = start; i < end; i++)
|
||||
size_t start = patch ? patch-1 : 0;
|
||||
size_t end = patch ? start : myModel.size();
|
||||
for (size_t i = start; i < end && i < myModel.size(); i++)
|
||||
if (!myModel[i]->empty())
|
||||
if (!myModel[i]->write(os,basis))
|
||||
return false;
|
||||
@ -1438,7 +1513,7 @@ bool SIMbase::dumpSolution (const Vector& psol, std::ostream& os) const
|
||||
|
||||
|
||||
bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
|
||||
bool formatted, std::streamsize outputPrecision) const
|
||||
bool formatted, std::streamsize precision) const
|
||||
{
|
||||
if (psol.empty() || myPoints.empty())
|
||||
return true;
|
||||
@ -1491,8 +1566,8 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
|
||||
return false;
|
||||
|
||||
// Formatted output, use scientific notation with fixed field width
|
||||
std::streamsize flWidth = 8 + outputPrecision;
|
||||
std::streamsize oldPrec = os.precision(outputPrecision);
|
||||
std::streamsize flWidth = 8 + precision;
|
||||
std::streamsize oldPrec = os.precision(precision);
|
||||
std::ios::fmtflags oldF = os.flags(std::ios::scientific | std::ios::right);
|
||||
for (j = 0; j < points.size(); j++)
|
||||
{
|
||||
@ -1537,21 +1612,81 @@ bool SIMbase::dumpResults (const Vector& psol, double time, std::ostream& os,
|
||||
}
|
||||
|
||||
|
||||
bool SIMbase::project (Vector& psol)
|
||||
bool SIMbase::project (Matrix& ssol, const Vector& psol,
|
||||
ProjectionMethod pMethod) const
|
||||
{
|
||||
Matrix values;
|
||||
Vector ssol;
|
||||
ssol.resize(psol.size()*myProblem->getNoFields());
|
||||
for (size_t i = 0; i < myModel.size(); i++)
|
||||
{
|
||||
myModel[i]->extractNodeVec(psol,myProblem->getSolution());
|
||||
if (!myModel[i]->evalSolution(values,*myProblem))
|
||||
return false;
|
||||
if (!myModel[i]->injectNodeVec(values,ssol,values.rows()))
|
||||
return false;
|
||||
}
|
||||
psol = ssol;
|
||||
PROFILE1("Solution projection");
|
||||
|
||||
if (msgLevel > 1)
|
||||
std::cout <<"\nProjecting secondary solution ...\n"<< std::endl;
|
||||
|
||||
ssol.resize(0,0);
|
||||
|
||||
size_t i, j, n;
|
||||
size_t ngNodes = mySam->getNoNodes();
|
||||
|
||||
Matrix values;
|
||||
Vector count(myModel.size() > 1 ? ngNodes : 0);
|
||||
|
||||
for (i = 0; i < myModel.size(); i++)
|
||||
{
|
||||
if (myModel[i]->empty()) continue; // skip empty patches
|
||||
|
||||
// Extract the primary solution control point values for this patch
|
||||
myModel[i]->extractNodeVec(psol,myProblem->getSolution());
|
||||
|
||||
// Project the secondary solution and retrieve control point values
|
||||
switch (pMethod) {
|
||||
case GLOBAL:
|
||||
if (!myModel[i]->evalSolution(values,*myProblem))
|
||||
return false;
|
||||
break;
|
||||
|
||||
case LOCAL:
|
||||
// Annette, add your local projection stuff here...
|
||||
|
||||
default:
|
||||
std::cerr <<" *** SIMbase::project: Projection method "<< pMethod
|
||||
<<" not implemented."<< std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
size_t nComps = values.rows();
|
||||
size_t nNodes = myModel[i]->getNoNodes();
|
||||
if (ssol.empty())
|
||||
ssol.resize(nComps,ngNodes);
|
||||
|
||||
// Nodal averaging for nodes referred to by two or more patches
|
||||
// (these are typically the interface nodes)
|
||||
for (n = 1; n <= nNodes; n++)
|
||||
if (count.empty())
|
||||
ssol.fillColumn(myModel[i]->getNodeID(n),values.getColumn(n));
|
||||
else
|
||||
{
|
||||
int inod = myModel[i]->getNodeID(n);
|
||||
for (j = 1; j <= nComps; j++)
|
||||
ssol(j,inod) += values(j,n);
|
||||
count(inod) ++;
|
||||
}
|
||||
}
|
||||
|
||||
// Divide through by count(n) to get the nodal average at the interface nodes
|
||||
for (n = 1; n <= count.size(); n++)
|
||||
if (count(n) > 1.0)
|
||||
for (j = 1; j <= ssol.rows(); j++)
|
||||
ssol(j,n) /= count(n);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool SIMbase::project (Vector& sol) const
|
||||
{
|
||||
Matrix secsol;
|
||||
if (!this->project(secsol,sol))
|
||||
return false;
|
||||
|
||||
sol = secsol;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
@ -90,16 +90,16 @@ public:
|
||||
//! \param is The file stream to read from
|
||||
virtual bool parse(char* keyWord, std::istream& is);
|
||||
|
||||
//! \brief Reads patches from a given stream.
|
||||
//! \param[in] isp The stream to read from
|
||||
virtual void readPatches(std::istream& isp) {}
|
||||
|
||||
//! \brief Performs some pre-processing tasks on the FE model.
|
||||
//! \param[in] ignoredPatches Indices of patches to ignore in the analysis
|
||||
//! \param[in] fixDup Merge duplicated FE nodes on patch interfaces?
|
||||
bool preprocess(const std::vector<int>& ignoredPatches = std::vector<int>(),
|
||||
bool fixDup = false);
|
||||
|
||||
//! \brief Read patches from given stream
|
||||
//! \param[in] isp The stream to read from
|
||||
virtual void readPatches(std::istream& isp) {};
|
||||
|
||||
//! \brief Allocates the system matrices of the FE problem to be solved.
|
||||
//! \param[in] mType The matrix format to use
|
||||
//! \param[in] nMats Number of system matrices
|
||||
@ -215,24 +215,47 @@ public:
|
||||
|
||||
//! \brief Integrates some solution norm quantities.
|
||||
//! \details If an analytical solution is provided, norms of the exact
|
||||
//! error in the solution are computed as well.
|
||||
//! error in the solution are computed as well. If projected secondary
|
||||
//! solutions are provided (i.e., \a ssol is not empty), norms of the
|
||||
//! difference between these solutions and the directly evaluated secondary
|
||||
//! solution are computed as well.
|
||||
//! \param[in] time Parameters for nonlinear/time-dependent simulations.
|
||||
//! \param[in] psol Primary solution vectors
|
||||
//! \param[in] ssol Secondary solution vectors
|
||||
//! \param[out] gNorm Global norm quantities
|
||||
//! \param[out] eNorm Element-wise norm quantities
|
||||
bool solutionNorms(const TimeDomain& time,
|
||||
const Vectors& psol, const Vectors& ssol,
|
||||
Vector& gNorm, Matrix* eNorm = 0);
|
||||
//! \brief Integrates some solution norm quantities.
|
||||
//! \param[in] time Parameters for nonlinear/time-dependent simulations.
|
||||
//! \param[in] psol Primary solution vectors
|
||||
//! \param[out] gNorm Global norm quantities
|
||||
//! \param[out] eNorm Element-wise norm quantities
|
||||
virtual bool solutionNorms(const TimeDomain& time, const Vectors& psol,
|
||||
Vector& gNorm, Matrix* eNorm = 0);
|
||||
|
||||
//!
|
||||
//! \details Use this version if no projected solutions are needed/available.
|
||||
bool solutionNorms(const TimeDomain& time, const Vectors& psol,
|
||||
Vector& gNorm, Matrix* eNorm = 0)
|
||||
{ return this->solutionNorms(time,psol,Vectors(),gNorm,eNorm); }
|
||||
//! \brief Integrates some solution norm quantities.
|
||||
//! \details If an analytical solution is provided, norms of the exact
|
||||
//! error in the solution are computed as well.
|
||||
//! \param[in] psol Primary solution vectors
|
||||
//! \param[in] ssol Secondary solution vectors
|
||||
//! \param[out] eNorm Element-wise norm quantities
|
||||
//! \param[out] gNorm Global norm quantities
|
||||
//!
|
||||
//! \details Use this version for linear/stationary problems only.
|
||||
virtual bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
|
||||
{ return this->solutionNorms(TimeDomain(),psol,gNorm,&eNorm); }
|
||||
bool solutionNorms(const Vectors& psol, const Vectors& ssol,
|
||||
Matrix& eNorm, Vector& gNorm)
|
||||
{ return this->solutionNorms(TimeDomain(),psol,ssol,gNorm,&eNorm); }
|
||||
//! \brief Integrates some solution norm quantities.
|
||||
//! \param[in] psol Primary solution vectors
|
||||
//! \param[out] eNorm Element-wise norm quantities
|
||||
//! \param[out] gNorm Global norm quantities
|
||||
//!
|
||||
//! \details Use this version for linear/stationary problems,
|
||||
//! and when no projected solutions are needed/available.
|
||||
bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
|
||||
{ return this->solutionNorms(TimeDomain(),psol,Vectors(),gNorm,&eNorm); }
|
||||
|
||||
//! \brief Computes the total reaction forces in the model.
|
||||
//! \param[out] RF Reaction force in each spatial direction + energy
|
||||
@ -251,13 +274,21 @@ public:
|
||||
int nev, int ncv, int iop, double shift,
|
||||
size_t iA = 0, size_t iB = 1);
|
||||
|
||||
//! \brief Projects the secondary solution associated with \a psol.
|
||||
//! \param psol Control point values of primary(in)/secondary(out) solution
|
||||
//! \brief Enum defining the available projection methods.
|
||||
enum ProjectionMethod { GLOBAL, LOCAL };
|
||||
//! \brief Projects the secondary solution associated with a primary solution.
|
||||
//! \param[out] ssol Control point values of the secondary solution
|
||||
//! \param[in] psol Control point values of the primary solution
|
||||
//! \param[in] pMethod Projection method to use
|
||||
//!
|
||||
//! \details The secondary solution, defined through the Integrand object,
|
||||
//! corresponding to the primary solution \a psol is projected onto the
|
||||
//! spline basis to obtain the control point values of the secondary solution.
|
||||
bool project(Vector& psol);
|
||||
bool project(Matrix& ssol, const Vector& psol,
|
||||
ProjectionMethod pMethod = GLOBAL) const;
|
||||
//! \brief Projects the secondary solution associated with a primary solution.
|
||||
//! \param sol Control point values of primary(in)/secondary(out) solution
|
||||
bool project(Vector& sol) const;
|
||||
|
||||
//! \brief Evaluates the secondary solution field for specified patch.
|
||||
//! \param[out] field Control point values of the secondary solution field
|
||||
@ -324,6 +355,18 @@ public:
|
||||
const char* pvecName = 0, int idBlock = 10,
|
||||
int psolComps = 0);
|
||||
|
||||
//! \brief Writes projected solutions for a given time step to the VTF-file.
|
||||
//! \param[in] ssol Secondary solution vector (control point values)
|
||||
//! \param[in] nViz Number of visualization points over each knot-span
|
||||
//! \param[in] iStep Load/time step identifier
|
||||
//! \param nBlock Running result block counter
|
||||
//! \param[in] time Load/time step parameter
|
||||
//! \param[in] idBlock Starting value of result block numbering
|
||||
//! \param[in] prefix Common prefix for the field components
|
||||
bool writeGlvP(const Vector& ssol, const int* nViz, int iStep,
|
||||
int& nBlock, double time = 0.0, int idBlock = 100,
|
||||
const char* pvecName = "Global projected");
|
||||
|
||||
//! \brief Writes a mode shape and associated eigenvalue to the VTF-file.
|
||||
//! \details The eigenvalue is used as a label on the step state info
|
||||
//! that is associated with the eigenvector.
|
||||
@ -344,7 +387,7 @@ public:
|
||||
//! \brief Writes time/load step info to the VTF-file.
|
||||
//! \param[in] iStep Load/time step identifier
|
||||
//! \param[in] value Time or load parameter of the step
|
||||
//! \param[in] itype Type indentifier of the step
|
||||
//! \param[in] itype Type identifier of the step
|
||||
bool writeGlvStep(int iStep, double value = 0.0, int itype = 0);
|
||||
|
||||
//! \brief Closes the current VTF-file.
|
||||
@ -357,7 +400,7 @@ public:
|
||||
//! \param os Output stream to write the spline data to
|
||||
//! \param[in] basis Which basis to dump for mixed methods (0 = geometry)
|
||||
//! \param[in] patch Which patch to dump for (0 = all)
|
||||
bool dumpBasis(std::ostream& os, int basis = 0, size_t patch=0) const;
|
||||
bool dumpBasis(std::ostream& os, int basis = 0, size_t patch = 0) const;
|
||||
//! \brief Dumps the entire solution in ASCII format.
|
||||
//! \param[in] psol Primary solution vector to derive other quantities from
|
||||
//! \param os Output stream to write the solution data to
|
||||
@ -378,6 +421,9 @@ public:
|
||||
void dumpPrimSol(const Vector& psol, std::ostream& os,
|
||||
bool withID = true) const;
|
||||
|
||||
//! \brief Returns whether an analytical solution is available or not.
|
||||
bool haveAnaSol() const { return mySol ? true : false; }
|
||||
|
||||
protected:
|
||||
//! \brief Defines the type of a property set.
|
||||
//! \param[in] code The property code to be associated with the property type
|
||||
|
||||
Loading…
Reference in New Issue
Block a user