diff --git a/Apps/LinearElasticity/main_LinEl3D.C b/Apps/LinearElasticity/main_LinEl3D.C index 2aaaabe2..d49fda51 100644 --- a/Apps/LinearElasticity/main_LinEl3D.C +++ b/Apps/LinearElasticity/main_LinEl3D.C @@ -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|) : " + 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); diff --git a/src/ASM/IntegrandBase.C b/src/ASM/IntegrandBase.C index 80ab2f1c..ad265574 100644 --- a/src/ASM/IntegrandBase.C +++ b/src/ASM/IntegrandBase.C @@ -234,23 +234,37 @@ bool NormBase::initElementBou (const std::vector& 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(); } diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index 58346d9f..a4a02233 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -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; } diff --git a/src/Integrands/Elasticity.C b/src/Integrands/Elasticity.C index bb322f75..1a973d45 100644 --- a/src/Integrands/Elasticity.C +++ b/src/Integrands/Elasticity.C @@ -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; diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 734bcdbb..7165d468 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -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)) + { + const char* normName = NormBase::getName(j,this->haveAnaSol(),*p); + if (!myVtf->writeSblk(sID[k],normName,++idBlock,iStep,true)) return false; - else if (prefix && npc > 0) + + if (prefix && npc > 0) { - if (k == 2) + 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; } @@ -2155,13 +2159,14 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, if (!myModel[i]->evalSolution(values,*myProblem,NULL,'A')) return false; break; - + case QUASI: if (msgLevel > 1 && i == 0) std::cout <<"\tQuasi interpolation"<< std::endl; 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; diff --git a/src/Utility/HDF5Writer.C b/src/Utility/HDF5Writer.C index ad37f7c4..c13932e8 100644 --- a/src/Utility/HDF5Writer.C +++ b/src/Utility/HDF5Writer.C @@ -352,7 +352,7 @@ void HDF5Writer::writeSIM (int level, const DataEntry& entry, for (size_t j=0;jhaveAnaSol()),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;jhaveAnaSol()),0,&dummy,H5T_NATIVE_DOUBLE); } } } diff --git a/src/Utility/XMLWriter.C b/src/Utility/XMLWriter.C index 419fb1f3..4701e2b7 100644 --- a/src/Utility/XMLWriter.C +++ b/src/Utility/XMLWriter.C @@ -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"); - } + NormBase* norm = sim->getNormIntegrand(); + 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; } } @@ -208,5 +212,5 @@ int XMLWriter::realTimeLevel(int filelevel) const int XMLWriter::realTimeLevel(int filelevel, int order, int interval) const { - return filelevel/order*interval; + return filelevel/order*interval; }