From aeae2faad1ef21f59c8e2d19aadf1fb019695d63 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 26 Oct 2017 07:22:37 +0200 Subject: [PATCH] Changed: Use a lambda function in assembleSystem to reduce code duplication. Changed: solveMatrixSystem is renamed to solveSystem (overloaded method). Added: Support for dual solution field, for goal-oriented error estimates. --- src/SIM/AdaptiveSIM.C | 28 ++++++++++---- src/SIM/AdaptiveSIM.h | 3 ++ src/SIM/AdaptiveSetup.C | 10 ++++- src/SIM/AdaptiveSetup.h | 3 +- src/SIM/SIMbase.C | 82 +++++++++++++++++++++++------------------ src/SIM/SIMbase.h | 8 ++-- 6 files changed, 85 insertions(+), 49 deletions(-) diff --git a/src/SIM/AdaptiveSIM.C b/src/SIM/AdaptiveSIM.C index e431aeaf..5ae06b60 100644 --- a/src/SIM/AdaptiveSIM.C +++ b/src/SIM/AdaptiveSIM.C @@ -44,6 +44,9 @@ bool AdaptiveSIM::initAdaptor (size_t normGroup) return false; projs.resize(opt.project.size()); + if (model.haveDualSol()) + projd.resize(opt.project.size()); + if (opt.format >= 0) { prefix.reserve(opt.project.size()); @@ -101,7 +104,9 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF, }; gNorm.clear(); + dNorm.clear(); eNorm.clear(); + fNorm.clear(); model.getProcessAdm().cout <<"\nAdaptive step "<< iStep << std::endl; if (iStep > 1) @@ -127,17 +132,17 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF, if (!this->assembleAndSolveSystem()) return failure(); - eNorm.clear(); - gNorm.clear(); - // Project the secondary solution onto the splines basis size_t idx = 0; model.setMode(SIM::RECOVERY); for (const SIMoptions::ProjectionMap::value_type& prj : opt.project) if (prj.first <= SIMoptions::NONE) - projs[idx++].clear(); // No projection for this norm group + idx++; // No projection for this norm group else if (!model.project(projs[idx++],solution.front(),prj.first)) return failure(); + else if (idx == adaptor && idx <= projd.size() && solution.size() > 1) + if (!model.project(projd[idx-1],solution[1],prj.first)) + return failure(); if (msgLevel > 1 && !projs.empty()) model.getProcessAdm().cout << std::endl; @@ -147,6 +152,9 @@ bool AdaptiveSIM::solveStep (const char* inputfile, int iStep, bool withRF, model.setQuadratureRule(opt.nGauss[1]); if (!model.solutionNorms(solution.front(),projs,eNorm,gNorm)) return failure(); + if (!projd.empty() && solution.size() > 1) + if (!model.solutionNorms(solution[1],projd,fNorm,dNorm)) + return failure(); model.setMode(SIM::RECOVERY); if (!model.dumpResults(solution.front(),0.0, @@ -172,15 +180,21 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec) if (outPrec > 0) { std::streamsize oldPrec = IFEM::cout.precision(outPrec); - this->printNorms(gNorm,eNorm); + this->printNorms(gNorm,dNorm,eNorm); IFEM::cout.precision(oldPrec); } else - this->printNorms(gNorm,eNorm); + this->printNorms(gNorm,dNorm,eNorm); + + // Calculate refinement indicators including the dual error estimates + Vector refIn = eNorm.getRow(eRow); + if (eRow <= fNorm.rows()) + for (size_t i = 1; i <= refIn.size(); i++) + refIn(i) = sqrt(refIn(i)*fNorm(eRow,i)); // Set up refinement parameters LR::RefineData prm; - if (this->calcRefinement(prm,iStep,gNorm,eNorm.getRow(eRow)) <= 0) + if (this->calcRefinement(prm,iStep,gNorm,refIn) <= 0) return false; // Now refine the mesh and write out resulting grid diff --git a/src/SIM/AdaptiveSIM.h b/src/SIM/AdaptiveSIM.h index 698d130e..c91124ff 100644 --- a/src/SIM/AdaptiveSIM.h +++ b/src/SIM/AdaptiveSIM.h @@ -79,12 +79,15 @@ protected: private: Vectors gNorm; //!< Global norms + Vectors dNorm; //!< Dual global norms Matrix eNorm; //!< Element norms + Matrix fNorm; //!< Dual element norms int geoBlk; //!< Running VTF geometry block counter int nBlock; //!< Running VTF result block counter std::vector projs; //!< Projected secondary solutions + std::vector projd; //!< Projected dual solutions std::vector prefix; //!< Norm prefices for VTF-output protected: diff --git a/src/SIM/AdaptiveSetup.C b/src/SIM/AdaptiveSetup.C index 22d044d2..6baec3f2 100644 --- a/src/SIM/AdaptiveSetup.C +++ b/src/SIM/AdaptiveSetup.C @@ -438,7 +438,7 @@ int AdaptiveSetup::calcRefinement (LR::RefineData& prm, int iStep, } -void AdaptiveSetup::printNorms (const Vectors& gNorm, +void AdaptiveSetup::printNorms (const Vectors& gNorm, const Vectors& dNorm, const Matrix& eNorm, size_t w) const { model.printNorms(gNorm,w); @@ -457,6 +457,14 @@ void AdaptiveSetup::printNorms (const Vectors& gNorm, if (model.haveAnaSol() && adaptor != 0) IFEM::cout <<"\nEffectivity index : " << model.getEffectivityIndex(gNorm,adaptor,adNorm); + if (adaptor < dNorm.size()) + { + const Vector& bNorm = dNorm[adaptor]; + IFEM::cout <<"\nError estimate, dual solution" + << utl::adjustRight(w-29,"(z)") << bNorm(adNorm) + <<"\nRelative error (%) : " + << 100.0*bNorm(adNorm)/hypot(dNorm.front()(1),bNorm(2)); + } IFEM::cout <<"\n"<< std::endl; delete norm; diff --git a/src/SIM/AdaptiveSetup.h b/src/SIM/AdaptiveSetup.h index 49fed1f6..a35ea95f 100644 --- a/src/SIM/AdaptiveSetup.h +++ b/src/SIM/AdaptiveSetup.h @@ -69,9 +69,10 @@ public: //! \brief Prints out global norms to the log stream. //! \param[in] gNorm Global norms + //! \param[in] dNorm Global dual norms //! \param[in] eNorm Element norms //! \param[in] w Field width for global norm labels - void printNorms(const Vectors& gNorm, + void printNorms(const Vectors& gNorm, const Vectors& dNorm, const Matrix& eNorm, size_t w = 36) const; //! \brief Dumps current mesh to external file(s) for inspection. diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 04bc1b03..eb913a85 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -661,7 +661,7 @@ size_t SIMbase::getNoEquations () const size_t SIMbase::getNoRHS () const { - return myEqSys ? myEqSys->getNoRHS() : 1; + return myEqSys ? myEqSys->getNoRHS() : (this->haveDualSol() ? 2 : 1); } @@ -809,6 +809,37 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol, { PROFILE1("Element assembly"); + // Lambda function for assembling the interior terms for a given patch + auto&& assembleInterior = [this,time,prevSol](IntegrandBase* integrand, + GlobalIntegral& integral, + ASMbase* pch, int pidx) + { + if (msgLevel > 1) + IFEM::cout <<"\nAssembling interior matrix terms for P"<< pidx + << std::endl; + + if (!this->initBodyLoad(pidx)) + return false; + + if (!this->extractPatchSolution(integrand,prevSol,pidx-1)) + return false; + + if (dualField && myProblem->getExtractionField()) + { + Matrix extrField(*myProblem->getExtractionField()); + if (!pch->L2projection(extrField,dualField)) + return false; + } + + if (mySol) + mySol->initPatch(pch->idx); + +#if SP_DEBUG > 2 + integrand->printSolution(std::cout,pch->idx+1); +#endif + return pch->integrate(*integrand,integral,time); + }; + bool ok = true; bool isAssembling = (myProblem->getMode() > SIM::INIT && myProblem->getMode() < SIM::RECOVERY); @@ -840,6 +871,7 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol, for (p = myProps.begin(); p != myProps.end() && ok; ++p) if (p->pcode == Property::MATERIAL && (it->first == 0 || it->first == p->pindx)) + { if (!(pch = this->getPatch(p->patch))) { std::cerr <<" *** SIMbase::assembleSystem: Patch index "<< p->patch @@ -849,39 +881,17 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol, else if (this->initMaterial(p->pindx)) { lp = p->patch; - if (msgLevel > 1) - IFEM::cout <<"\nAssembling interior matrix terms for P"<< lp - << std::endl; - ok &= this->initBodyLoad(lp); - ok &= this->extractPatchSolution(it->second,prevSol,lp-1); - if (mySol) - mySol->initPatch(pch->idx); -#if SP_DEBUG > 2 - it->second->printSolution(std::cout,pch->idx+1); -#endif - ok &= pch->integrate(*it->second,sysQ,time); + ok = assembleInterior(it->second,sysQ,pch,lp); } else ok = false; + } if (lp == 0 && it->first == 0) // All patches refer to the same material, and we assume it has been // initialized during input processing (thus no initMaterial call here) - for (size_t k = 0; k < myModel.size() && ok; k++) - { - lp = k+1; - if (msgLevel > 1) - IFEM::cout <<"\nAssembling interior matrix terms for P"<< lp - << std::endl; - ok &= this->initBodyLoad(lp); - ok &= this->extractPatchSolution(it->second,prevSol,k); - if (mySol) - mySol->initPatch(myModel[k]->idx); -#if SP_DEBUG > 2 - it->second->printSolution(std::cout,myModel[k]->idx+1); -#endif - ok &= myModel[k]->integrate(*it->second,sysQ,time); - } + for (lp = 0; lp < myModel.size() && ok; lp++) + ok = assembleInterior(it->second,sysQ,myModel[lp],1+lp); } // Assemble contributions from the Neumann boundary conditions @@ -1087,17 +1097,17 @@ bool SIMbase::solveSystem (Vector& solution, int printSol, double* rCond, } -bool SIMbase::solveMatrixSystem (Vectors& solution, int printSol, - const char* compName) +bool SIMbase::solveSystem (Vectors& solution, int printSol, const char* cmpName) { - solution.resize(this->getNoRHS()); - for (size_t i = 0; i < solution.size(); i++) - if (!this->solveSystem(solution[i],printSol,nullptr,compName,i==0,i)) - return false; - else if (solution.size() > 2) - printSol = 0; // Print summary only for the first two solutions + size_t nSol = myEqSys ? myEqSys->getNoRHS() : 0; + if (solution.size() < nSol) + solution.resize(nSol); - return true; + bool status = nSol > 0; + for (size_t i = 0; i < nSol && status; i++) + status = this->solveSystem(solution[i],printSol,nullptr,cmpName,i==0,i); + + return status; } diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index 8a2988b8..dc1e9218 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -285,11 +285,11 @@ public: } //! \brief Solves a linear system of equations with multiple right-hand-sides. - //! \param[out] solution Global primary solution vector + //! \param[out] solution Global primary solution vectors //! \param[in] printSol Print solution if its size is less than \a printSol - //! \param[in] compName Solution name to be used in norm output - bool solveMatrixSystem(Vectors& solution, int printSol = 0, - const char* compName = "displacement"); + //! \param[in] cmpName Solution name to be used in norm output + bool solveSystem(Vectors& solution, int printSol = 0, + const char* cmpName = "displacement"); //! \brief Finds the DOFs showing the worst convergence behavior. //! \param[in] x Global primary solution vector