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

Changed: Move the calculation of the total number of integration points and the associated buffer allocation until after the quadrature rule has been defined

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1426 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2012-01-22 18:12:45 +00:00
committed by Knut Morten Okstad
parent e5b51b0b26
commit 25de55cc65
5 changed files with 52 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/SIM/SIM1D.C
View File

@@ -17,7 +17,6 @@
#include "ASMs1DSpec.h"
#include "Functions.h"
#include "Utilities.h"
#include "tinyxml.h"
@@ -411,6 +410,8 @@ void SIM1D::setQuadratureRule (size_t ng)
for (size_t i = 0; i < myModel.size(); i++)
if (!myModel.empty())
static_cast<ASMs1D*>(myModel[i])->setGauss(ng);
this->initIntegrationBuffers();
}

10
src/SIM/SIM2D.C
View File

@@ -15,9 +15,7 @@
#include "ASMs2DC1.h"
#include "Functions.h"
#include "Utilities.h"
#include "tinyxml.h"
#ifdef USE_OPENMP
#include <omp.h>
#endif
@@ -203,7 +201,7 @@ bool SIM2D::parseGeometryTag(const TiXmlElement* elem)
std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs2D*>(myModel[patch-1])->closeEdges(pedir);
// cannot do multi-threaded assembly with periodicities
// Cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
@@ -390,11 +388,11 @@ bool SIM2D::parse (char* keyWord, std::istream& is)
std::cout <<"\tPeriodic "<< char('H'+pedir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs2D*>(myModel[patch-1])->closeEdges(pedir);
// cannot do multi-threaded assembly with periodicities
}
// Cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
}
}
else if (!strncasecmp(keyWord,"CONSTRAINTS",11))
@@ -552,6 +550,8 @@ void SIM2D::setQuadratureRule (size_t ng)
for (size_t i = 0; i < myModel.size(); i++)
if (!myModel.empty())
static_cast<ASMs2D*>(myModel[i])->setGauss(ng);
this->initIntegrationBuffers();
}

13
src/SIM/SIM3D.C
View File

@@ -17,9 +17,11 @@
#include "ASMs3DSpec.h"
#include "Functions.h"
#include "Utilities.h"
#include <fstream>
#include "tinyxml.h"
#include <fstream>
#ifdef USE_OPENMP
#include <omp.h>
#endif
SIM3D::SIM3D (bool checkRHS, unsigned char n1, unsigned char n2)
@@ -211,11 +213,12 @@ bool SIM3D::parse (char* keyWord, std::istream& is)
std::cout <<"\tPeriodic "<< char('H'+pfdir) <<"-direction P"<< patch
<< std::endl;
static_cast<ASMs3D*>(myModel[patch-1])->closeFaces(pfdir);
// cannot do multi-threaded assembly with periodicities
}
// Cannot do multi-threaded assembly with periodicities
#ifdef USE_OPENMP
omp_set_num_threads(1);
#endif
}
}
else if (!strncasecmp(keyWord,"CONSTRAINTS",11))
@@ -565,6 +568,8 @@ void SIM3D::setQuadratureRule (size_t ng)
for (size_t i = 0; i < myModel.size(); i++)
if (!myModel.empty())
static_cast<ASMs3D*>(myModel[i])->setGauss(ng);
this->initIntegrationBuffers();
}

49
src/SIM/SIMbase.C
View File

@@ -58,6 +58,7 @@ SIMbase::SIMbase () : g2l(&myGlb2Loc)
vizIncr = 1;
format = 1;
mixedFEM = false;
nIntGP = nBouGP = 0;
MPCLess::compareSlaveDofOnly = true; // to avoid multiple slave definitions
@@ -658,13 +659,9 @@ bool SIMbase::preprocess (const std::vector<int>& ignored, bool fixDup)
if (renum > 0)
std::cout <<"\nRenumbered "<< renum <<" nodes."<< std::endl;
size_t nIntP = 0, nBouP = 0;
ASMs2DC1::renumberNodes(*g2l);
for (mit = myModel.begin(); mit != myModel.end(); mit++)
{
(*mit)->renumberNodes(*g2l);
(*mit)->getNoIntPoints(nIntP); // count number of global integration points
}
// Process the specified Dirichlet boundary conditions
std::cout <<"\nResolving Dirichlet boundary conditions"<< std::endl;
@@ -687,12 +684,6 @@ bool SIMbase::preprocess (const std::vector<int>& ignored, bool fixDup)
ok = false;
break;
case Property::NEUMANN:
// Count number of global boundary integration points
if (p->patch > 0 && p->patch <= myModel.size())
myModel[p->patch-1]->getNoBouPoints(nBouP,p->ldim,p->lindx);
break;
default:
break;
}
@@ -709,9 +700,6 @@ bool SIMbase::preprocess (const std::vector<int>& ignored, bool fixDup)
// Resolve possibly chaining of the MPC equations
if (!allMPCs.empty()) ASMbase::resolveMPCchains(allMPCs);
// Let the integrand know how many integration points in total do we have
if (myProblem) myProblem->initIntegration(nIntP,nBouP);
// Preprocess the result points
for (ResPointVec::iterator p = myPoints.begin(); p != myPoints.end();)
{
@@ -753,6 +741,23 @@ bool SIMbase::preprocess (const std::vector<int>& ignored, bool fixDup)
}
void SIMbase::initIntegrationBuffers ()
{
nIntGP = nBouGP = 0;
for (size_t i = 0; i < myModel.size(); i++)
myModel[i]->getNoIntPoints(nIntGP); // count the interior integration points
for (PropertyVec::const_iterator p = myProps.begin(); p != myProps.end(); p++)
if (p->pcode == Property::NEUMANN)
// Count the boundary integration points
if (p->patch > 0 && p->patch <= myModel.size())
myModel[p->patch-1]->getNoBouPoints(nBouGP,p->ldim,p->lindx);
// Let the integrand know how many integration points in total do we have
if (myProblem) myProblem->initIntegration(nIntGP,nBouGP);
}
bool SIMbase::setPropertyType (int code, Property::Type ptype, int pindex)
{
if (code < 0)
@@ -762,12 +767,11 @@ bool SIMbase::setPropertyType (int code, Property::Type ptype, int pindex)
return false;
}
for (size_t j = 0; j < myProps.size(); j++)
if (myProps[j].pindx == (size_t)code &&
myProps[j].pcode == Property::UNDEFINED)
for (PropertyVec::iterator p = myProps.begin(); p != myProps.end(); p++)
if (p->pindx == (size_t)code && p->pcode == Property::UNDEFINED)
{
myProps[j].pcode = ptype;
if (pindex >= 0) myProps[j].pindx = pindex;
p->pcode = ptype;
if (pindex >= 0) p->pindx = pindex;
}
return true;
@@ -948,7 +952,7 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
else
ok = false;
if (j == 0 && ok)
if (j == 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 (i = 0; i < myModel.size() && ok; i++)
@@ -1001,6 +1005,8 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
else
ok = false;
if (!ok) std::cerr <<" *** SIMbase::assembleSystem: Failure.\n"<< std::endl;
return ok && this->finalizeAssembly(newLHSmatrix);
}
@@ -1162,6 +1168,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
}
myProblem->initIntegration(time);
norm->initIntegration(nIntGP,nBouGP);
const Vector& primsol = psol.front();
size_t i, j, k;
@@ -1214,7 +1221,7 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
else
ok = false;
if (j == 0 && ok)
if (j == 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 (i = 0; i < myModel.size() && ok; i++)
@@ -1257,6 +1264,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
else
ok = false;
if (!ok) std::cerr <<" *** SIMbase::solutionNorms: Failure.\n"<< std::endl;
// Clean up the dynamically allocated norm objects. This will also perform
// the actual global norm assembly, in case the element norms are stored,
// and always when multi-threading is used.

7
src/SIM/SIMbase.h
View File

@@ -470,6 +470,9 @@ public:
NormBase* getNormIntegrand() const;
protected:
//! \brief Allocates the problem-dependent integration point buffers, if any.
void initIntegrationBuffers();
//! \brief Defines the type of a property set.
//! \param[in] code The property code to be associated with the property type
//! \param[in] ptype The property type to be associated with the given code
@@ -614,7 +617,11 @@ protected:
AlgEqSystem* myEqSys; //!< The actual linear equation system
SAMpatch* mySam; //!< Auxiliary data for FE assembly management
LinSolParams* mySolParams; //!< Input parameters for PETSc
private:
size_t nIntGP; //!< Number of interior integration points in the whole model
size_t nBouGP; //!< Number of boundary integration points in the whole model
//! \brief Parses a subelement of the <geometry> tag from input file
//! \param[in] elem The XML element to parse
bool parseGeometryTag(const TiXmlElement* elem);