OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixed: Corrected the norm output to VTF. It has been wrong since r1501.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1505 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2012-03-06 19:00:42 +00:00 committed by Knut Morten Okstad
parent e4bc621c31
commit aa5c798988
7 changed files with 65 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Apps/LinearElasticity/main_LinEl3D.C
View File

@ -68,9 +68,9 @@
\arg -DGL2 : Estimate error using discrete global L2 projection
\arg -CGL2 : Estimate error using continuous global L2 projection
\arg -SCR : Estimate error using Superconvergent recovery at Greville points
\arg -VDSA: Generate visualization fields through Variational Diminishing Spline Approximations
\arg -LSQ : Generate visualization fields through Least Square projections
\arg -QUASI : Generate visualization fields through Quasi interpolation projections
\arg -VDSA: Estimate error using Variational Diminishing Spline Approximations
\arg -LSQ : Estimate error using through Least Square projections
\arg -QUASI : Estimate error using Quasi-interpolation projections
*/
int main (int argc, char** argv)
@ -342,23 +342,24 @@ int main (int argc, char** argv)
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;
<<"\nExact relative error (%) : "<< gNorm(4)/gNorm(3)*100.0
<< std::endl;
size_t j = model->haveAnaSol() ? 5 : 3;
for (pit = pOpt.begin(); pit != pOpt.end() && j < gNorm.size(); pit++)
{
std::cout <<"\n>>> Error estimates based on "<< pit->second <<" <<<";
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;
std::cout <<"\nL2-norm, |u^r| : "<< gNorm(j++);
std::cout <<"\nL2-Error norm, 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 <<"\nL2-norm |s^r| =(s^r,s^r)^0.5 : "<< gNorm(j++);
std::cout <<"\nL2-error (e,e)^0.5, e=s^r-s^h : "<< gNorm(j++);
std::cout <<"\n- relative error (% of |s^r|) : "
<< gNorm(j-1)/gNorm(j-2)*100.0;
std::cout << std::endl;
}
@ -470,8 +471,7 @@ int main (int argc, char** argv)
return 12;
// Write element norms
if (!model->writeGlvN(eNorm,iStep,nBlock,prefix,
model->haveAnaSol() ? 3 : 2))
if (!model->writeGlvN(eNorm,iStep,nBlock,prefix,5))
return 13;
model->writeGlvStep(1);

28
src/ASM/IntegrandBase.C
View File

@ -234,23 +234,37 @@ bool NormBase::initElementBou (const std::vector<int>& MNPC1,
}
const char* NormBase::getName (size_t i, const char* prefix)
const char* NormBase::getName (size_t& i, bool haveExact, const char* prefix)
{
static const char* s[9] = {
static const char* s[8] = {
"a(u^h,u^h)^0.5",
"a(u,u)^0.5",
"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"
"L2(s^r)",
"L2(e), e=s^r-s^h"
};
if (!prefix) return s[i];
if (!haveExact) // Skip norms present only together with exact solutions
if (i == 1)
i = 3;
else if (i == 5)
i = 6;
if (!prefix && i < 8) return s[i];
static std::string name;
name = prefix + std::string(" ") + s[i];
if (prefix) name = prefix + std::string(" ");
if (i < 8)
name += s[i];
else
{
char comp[16];
sprintf(comp,"norm_%lu",i);
name += comp;
}
return name.c_str();
}

2
src/ASM/IntegrandBase.h
View File

@ -285,7 +285,7 @@ public:
virtual size_t getNoFields() const { return 0; }
//! \brief Returns the name of a norm quantity.
static const char* getName(size_t, const char* = 0);
static const char* getName(size_t&, bool, const char* = 0);
//! \brief Returns whether a norm component stores element contributions.
static bool hasElementContributions(size_t i) { return i != 1; }

17
src/Integrands/Elasticity.C
View File

@ -643,10 +643,10 @@ ElasticityNorm::ElasticityNorm (Elasticity& p, STensorFunc* a)
size_t ElasticityNorm::getNoFields () const
{
size_t nf = anasol ? 6 : 2;
size_t nf = anasol ? 4 : 2;
for (size_t i = 0; i < prjsol.size(); i++)
if (!prjsol.empty())
nf += anasol ? 5 : 2;
nf += anasol ? 5 : 4;
return nf;
}
@ -718,12 +718,8 @@ bool ElasticityNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
error = sigmar - sigmah;
pnorm[ip++] += error.dot(Cinv*error)*detJW;
// Integrate the L2-norm ||S|| = Int_Omega S:S dV
double tmp2 = sigmar.norm2();
pnorm[ip++] += tmp2*tmp2*detJW;
double tmp = error.norm2();
pnorm[ip++] += tmp*tmp*detJW;
double l2u = sigmar.norm2();
double l2e = error.norm2();
if (anasol)
{
@ -731,6 +727,11 @@ bool ElasticityNorm::evalInt (LocalIntegral& elmInt, const FiniteElement& fe,
error = sigma - sigmar;
pnorm[ip++] += error.dot(Cinv*error)*detJW;
}
// Integrate the L2-norm (sigma^r,sigma^r)
pnorm[ip++] += l2u*l2u*detJW;
// Integrate the error in L2-norm (sigma^r-sigma^h,sigma^r-sigma^h)
pnorm[ip++] += l2e*l2e*detJW;
}
return true;

13
src/SIM/SIMbase.C
View File

@ -1859,11 +1859,14 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock,
const char* z = 0;
const char** p = &z;
for (j = k = 0; k < maxN && !sID[k].empty(); j++, k++)
if (!myVtf->writeSblk(sID[k],NormBase::getName(j,*p),++idBlock,iStep,true))
return false;
else if (prefix && npc > 0)
{
if (k == 2)
const char* normName = NormBase::getName(j,this->haveAnaSol(),*p);
if (!myVtf->writeSblk(sID[k],normName,++idBlock,iStep,true))
return false;
if (prefix && npc > 0)
{
if (k == (this->haveAnaSol() ? 2 : 0))
p = prefix;
else if (j == 2+npc)
{
@ -1871,6 +1874,7 @@ bool SIMbase::writeGlvN (const Matrix& norms, int iStep, int& nBlock,
if (*p) p++;
}
}
}
return true;
}
@ -2162,6 +2166,7 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol,
if (!myModel[i]->evalSolution(values,*myProblem,NULL,'L'))
return false;
break;
case LEASTSQ:
if (msgLevel > 1 && i == 0)
std::cout <<"\tLeast squares projection"<< std::endl;

4
src/Utility/HDF5Writer.C
View File

@ -352,7 +352,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
for (size_t j=0;j<eNorm.rows();++j) {
if (NormBase::hasElementContributions(j)) {
Vector k;
writeArray(group2,NormBase::getName(j),patchEnorm.cols(),
writeArray(group2,NormBase::getName(j,sim->haveAnaSol()),patchEnorm.cols(),
patchEnorm.getRow(1+j).ptr(),H5T_NATIVE_DOUBLE);
}
}
@ -374,7 +374,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry,
if (entry.second.results & DataExporter::NORMS) {
for (size_t j=0;j<eNorm.rows();++j) {
if (NormBase::hasElementContributions(j))
writeArray(group2,NormBase::getName(j),0,&dummy,H5T_NATIVE_DOUBLE);
writeArray(group2,NormBase::getName(j,sim->haveAnaSol()),0,&dummy,H5T_NATIVE_DOUBLE);
}
}
}

16
src/Utility/XMLWriter.C
View File

@ -127,6 +127,7 @@ bool XMLWriter::readSIM (int level, const DataEntry& entry)
return true;
}
void XMLWriter::writeSIM (int level, const DataEntry& entry,
bool geometryUpdated)
{
@ -158,18 +159,21 @@ void XMLWriter::writeSIM (int level, const DataEntry& entry,
}
// secondary solution fields
size_t i, j;
if (entry.second.results & DataExporter::SECONDARY)
for (size_t j = 0; j < prob->getNoFields(2); j++)
addField(prob->getField2Name(j),"secondary",sim->getName()+(prob->mixedFormulation()?"-2":"-1"),1,sim->getNoPatches());
for (j = 0; j < prob->getNoFields(2); j++)
addField(prob->getField2Name(j),"secondary",
sim->getName() + (prob->mixedFormulation() ? "-2" : "-1"),
1,sim->getNoPatches());
// norms
if (entry.second.results & DataExporter::NORMS) {
// since the norm data isn't available, we have to instance the object
NormBase* norm = sim->getNormIntegrand();
for (size_t j = 0; j < norm->getNoFields(); j++) {
if (norm->hasElementContributions(j))
addField(norm->getName(j),"knotspan wise norm",sim->getName()+"-1",1,sim->getNoPatches(),"knotspan");
}
for (i = j = 0; i < norm->getNoFields(); i++, j++)
if (NormBase::hasElementContributions(i))
addField(NormBase::getName(j,sim->haveAnaSol()),"knotspan wise norm",
sim->getName()+"-1",1,sim->getNoPatches(),"knotspan");
delete norm;
}
}