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

Changed: Refactor globalL2projection into one common method that

Browse Source 
defines the projection matrices and solves the equation system, and a new
(patch-type dependent) method that carries out the assembly task.
Added: Support the continuous option in ASMs1D::globalL2projection.
This commit is contained in:
Knut Morten Okstad 2017-05-25 02:09:22 +02:00
parent 3ebdaca4d9
commit ad5bf4ed56
16 changed files with 256 additions and 265 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/ASM/ASMbase.C
View File

@ -1194,12 +1194,6 @@ bool ASMbase::evalSolution (Matrix&, const IntegrandBase&,
}
bool ASMbase::globalL2projection (Matrix&, const IntegrandBase&, bool) const
{
return Aerror("globalL2projection(Matrix&,const IntegrandBase&,bool)");
}
bool ASMbase::evaluate (const ASMbase*, const Vector&, RealArray&, int) const
{
return Aerror("evaluate(const ASMbase*,const Vector&,RealArray&,int)");

11
src/ASM/ASMbase.h
View File

@ -32,6 +32,8 @@ class GlobalIntegral;
class IntegrandBase;
class Integrand;
class ASMbase;
class SparseMatrix;
class StdVector;
class RealFunc;
class VecFunc;
class Vec3;
@ -623,6 +625,15 @@ protected:
int searchCtrlPt(RealArray::const_iterator cit, RealArray::const_iterator end,
const Vec3& X, int dimension, double tol = 0.001) const;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const = 0;
public:
//! \brief Constrains all nodes in the patch.
//! \param[in] dof Which DOFs to constrain at each node in the patch

79
src/ASM/ASMs1D.C
View File

@ -1359,42 +1359,34 @@ bool ASMs1D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
}
bool ASMs1D::globalL2projection (Matrix& sField,
const IntegrandBase& integrand,
bool continuous) const
bool ASMs1D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
if (!curv) return true; // silently ignore empty patches
if (continuous)
{
std::cerr <<" *** ASMs1D::globalL2projection: Only discrete L2-projection"
<<" is available in this method.\n "
<<" Use ASMbase::L2projection instead."<< std::endl;
}
const size_t nnod = this->getNoNodes();
const int p1 = curv->order();
// Get Gaussian quadrature point coordinates
const int ng1 = p1 - 1;
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr;
if (!xg) return false;
if (continuous && !wg) return false;
// Compute parameter values of the Gauss points over the whole patch
Matrix gp;
// and evaluate the secondary solution at all integration points
Matrix gp, sField;
RealArray gpar = this->getGaussPointParameters(gp,ng1,xg);
// Evaluate the secondary solution at all integration points
if (!this->evalSolution(sField,integrand,&gpar))
{
std::cerr <<" *** ASMs1D::assembleL2matrices: Failed for patch "<< idx+1
<<" nPoints="<< gpar.size() << std::endl;
return false;
}
// Set up the projection matrices
const size_t nnod = this->getNoNodes();
const size_t ncomp = sField.rows();
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod*ncomp);
A.redim(nnod,nnod);
double dL = 1.0;
Vector phi(p1);
Matrix dNdu, Xnod, J;
// === Assembly loop over all elements in the patch ==========================
@ -1404,12 +1396,32 @@ bool ASMs1D::globalL2projection (Matrix& sField,
{
if (MLGE[iel] < 1) continue; // zero-area element
if (continuous)
{
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,1+iel))
return false;
else if ((dL = this->getParametricLength(1+iel)) < 0.0)
return false; // topology error (probably logic error)
}
// --- Integration loop over all Gauss points in current element -----------
for (int i = 1; i <= ng1; i++, ip++)
for (int i = 0; i < ng1; i++, ip++)
{
// Fetch basis function values at current integration point
this->extractBasis(gp(i,1+iel),phi);
if (continuous)
this->extractBasis(gp(1+i,1+iel),phi,dNdu);
else
this->extractBasis(gp(1+i,1+iel),phi);
// Compute the Jacobian inverse and derivatives
double dJw = 1.0;
if (continuous)
{
dJw = 0.5*dL*wg[i]*utl::Jacobian(J,dNdu,Xnod,dNdu,false);
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the linear system A*x=B
for (size_t ii = 0; ii < phi.size(); ii++)
@ -1418,23 +1430,14 @@ bool ASMs1D::globalL2projection (Matrix& sField,
for (size_t jj = 0; jj < phi.size(); jj++)
{
int jnod = MNPC[iel][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj];
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= ncomp; r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1);
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
}
}
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
for (size_t i = 1; i <= nnod; i++)
for (size_t j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
return true;
}

20
src/ASM/ASMs1D.h
View File

@ -267,14 +267,6 @@ public:
//! \param[in] nSegSpan Number of visualization segments over each knot-span
virtual bool getGridParameters(RealArray& prm, int nSegSpan) const;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const;
//! \brief Establishes the vector with basis function values.
//! \param[in] u Parameter value of current integration point
//! \param[out] N Basis function values
@ -296,6 +288,15 @@ protected:
// Internal utility methods
// ========================
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const;
//! \brief Initializes the local element axes for a patch of beam elements.
//! \param[in] Zaxis Vector defining a point in the local XZ-plane
bool initLocalElementAxes(const Vec3& Zaxis);
@ -309,7 +310,8 @@ protected:
//! \param[in] master 0-based index of the first master node in this basis
//! \param[in] thick Thickness of connection
bool connectBasis(int vertex, ASMs1D& neighbor, int nvertex,
int basis = 1, int slave = 0, int master = 0, int thick = 1);
int basis = 1, int slave = 0, int master = 0,
int thick = 1);
//! \brief Extracts parameter values of the Gauss points.
//! \param[out] uGP Parameter values for all points

17
src/ASM/ASMs2D.h
View File

@ -461,19 +461,20 @@ public:
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool regular = true) const;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const;
protected:
// Internal utility methods
// ========================
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const;
//! \brief Connects all matching nodes on two adjacent boundary edges.
//! \param[in] edge Local edge index of this patch, in range [1,4]
//! \param neighbor The neighbor patch

40
src/ASM/ASMs2Drecovery.C
View File

@ -159,13 +159,11 @@ Go::GeomObject* ASMs2D::evalSolution (const IntegrandBase& integrand) const
}
bool ASMs2D::globalL2projection (Matrix& sField,
const IntegrandBase& integrand,
bool continuous) const
bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
if (!surf) return true; // silently ignore empty patches
PROFILE2("ASMs2D::globalL2");
const size_t nnod = this->getNoNodes(1);
const int p1 = surf->order_u();
const int p2 = surf->order_v();
@ -198,20 +196,14 @@ bool ASMs2D::globalL2projection (Matrix& sField,
surf->computeBasisGrid(gpar[0],gpar[1],spl0);
// Evaluate the secondary solution at all integration points
Matrix sField;
if (!this->evalSolution(sField,integrand,gpar.data()))
{
std::cerr <<" *** ASMs2D::globalL2projection: Failed for patch "<< idx+1
std::cerr <<" *** ASMs2D::assembleL2matrices: Failed for patch "<< idx+1
<<" nPoints="<< gpar[0].size()*gpar[1].size() << std::endl;
return false;
}
// Set up the projection matrices
const size_t nnod = this->getNoNodes(1);
const size_t ncomp = sField.rows();
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod*ncomp);
A.redim(nnod,nnod);
double dA = 1.0;
Vector phi(p1*p2);
Matrix dNdu, Xnod, J;
@ -262,30 +254,12 @@ bool ASMs2D::globalL2projection (Matrix& sField,
int jnod = MNPC[iel][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= ncomp; r++)
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
}
}
#if SP_DEBUG > 1
std::cout << " ---- Matrix A -----\n";
std::cout << A << std::endl;
std::cout << " -------------------\n";
std::cout << " ---- Vector B -----\n";
std::cout << B << std::endl;
std::cout << " -------------------\n";
#endif
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
for (size_t i = 1; i <= nnod; i++)
for (size_t j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
return true;
}

17
src/ASM/ASMs3D.h
View File

@ -524,19 +524,20 @@ public:
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool regular = true) const;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const;
protected:
// Internal utility methods
// ========================
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const;
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch

31
src/ASM/ASMs3Drecovery.C
View File

@ -159,13 +159,11 @@ Go::GeomObject* ASMs3D::evalSolution (const IntegrandBase& integrand) const
}
bool ASMs3D::globalL2projection (Matrix& sField,
const IntegrandBase& integrand,
bool continuous) const
bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
if (!svol) return true; // silently ignore empty patches
PROFILE2("ASMs3D::globalL2");
const size_t nnod = this->getNoNodes(1);
const int p1 = svol->order(0);
const int p2 = svol->order(1);
@ -204,21 +202,15 @@ bool ASMs3D::globalL2projection (Matrix& sField,
svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl0);
// Evaluate the secondary solution at all integration points
Matrix sField;
if (!this->evalSolution(sField,integrand,gpar.data()))
{
std::cerr <<" *** ASMs3D::globalL2projection: Failed for patch "<< idx+1
std::cerr <<" *** ASMs3D::assembleL2matrices: Failed for patch "<< idx+1
<<" nPoints="<< gpar[0].size()*gpar[1].size()*gpar[2].size()
<< std::endl;
return false;
}
// Set up the projection matrices
const size_t nnod = this->getNoNodes(1);
const size_t ncomp = sField.rows();
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod*ncomp);
A.redim(nnod,nnod);
double dV = 1.0;
Vector phi(p1*p2*p3);
Matrix dNdu, Xnod, J;
@ -271,21 +263,12 @@ bool ASMs3D::globalL2projection (Matrix& sField,
int jnod = MNPC[iel][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= ncomp; r++)
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
}
}
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
for (size_t i = 1; i <= nnod; i++)
for (size_t j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
return true;
}

80
src/ASM/GlbL2projector.C
View File

@ -37,11 +37,11 @@ public:
//! \brief Destruction method to clean up after numerical integration.
virtual void destruct() { delete elmData; delete this; }
GlbL2& gl2Int; //!< The global L2 projection integrand
LocalIntegral* elmData; //!< Element data associated with problem integrand
IntVec mnpc; //!< Matrix of element nodal correspondance for bases
std::vector<size_t> elem_sizes; //!< Size of each basis on the element
std::vector<size_t> basis_sizes; //!< Size of each basis on the patch
GlbL2& gl2Int; //!< The global L2-projection integrand
LocalIntegral* elmData; //!< Element data associated with problem integrand
IntVec mnpc; //!< Matrix of element nodal correspondance
std::vector<size_t> elem_sizes; //!< Size of each basis on the element
std::vector<size_t> basis_sizes; //!< Size of each basis on the patch
};
@ -63,7 +63,7 @@ LocalIntegral* GlbL2::getLocalIntegral (size_t nen, size_t iEl,
bool neumann) const
{
return new L2Mats(*const_cast<GlbL2*>(this),
problem.getLocalIntegral(nen,iEl,neumann));
problem.getLocalIntegral(nen,iEl,neumann));
}
@ -85,7 +85,7 @@ bool GlbL2::initElement (const IntVec& MNPC1,
{
L2Mats& gl2 = static_cast<L2Mats&>(elmInt);
gl2.mnpc = MNPC1;
gl2.mnpc = MNPC1;
gl2.elem_sizes = elem_sizes;
gl2.basis_sizes = basis_sizes;
return problem.initElement(MNPC1,elem_sizes,basis_sizes,*gl2.elmData);
@ -160,27 +160,35 @@ bool GlbL2::solve (Matrix& sField)
{
// Insert a 1.0 value on the diagonal for equations with no contributions.
// Needed in immersed boundary calculations with "totally outside" elements.
size_t nnod = A.dim();
for (size_t j = 1; j <= nnod; j++)
if (A(j,j) == 0.0) A(j,j) = 1.0;
size_t i, nnod = A.dim();
for (i = 1; i <= nnod; i++)
if (A(i,i) == 0.0) A(i,i) = 1.0;
#if SP_DEBUG > 1
std::cout <<"\nGlobal L2-projection matrix:\n"<< A;
std::cout <<"\nGlobal L2-projection RHS:"<< B;
#endif
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the nodal values of the projected field
size_t ncomp = B.dim() / nnod;
size_t j, ncomp = problem.getNoFields(2);
sField.resize(ncomp,nnod);
for (size_t i = 1; i <= nnod; i++)
for (size_t j = 1; j <= ncomp; j++)
for (i = 1; i <= nnod; i++)
for (j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
#if SP_DEBUG > 1
std::cout <<"\nSolution:"<< sField;
#endif
return true;
}
bool ASMbase::L2projection (Matrix& sField,
const IntegrandBase& integrand,
const TimeDomain& time)
const IntegrandBase& integrand,
const TimeDomain& time)
{
PROFILE2("ASMbase::L2projection");
@ -190,3 +198,45 @@ bool ASMbase::L2projection (Matrix& sField,
gl2.preAssemble(MNPC,this->getNoElms(true));
return this->integrate(gl2,dummy,time) && gl2.solve(sField);
}
bool ASMbase::globalL2projection (Matrix& sField,
const IntegrandBase& integrand,
bool continuous) const
{
if (this->empty()) return true; // silently ignore empty patches
PROFILE2("ASMbase::globalL2");
// Assemble the projection matrices
size_t i, nnod = this->getNoNodes(1);
size_t j, ncomp = integrand.getNoFields(2);
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod*ncomp);
A.redim(nnod,nnod);
if (!this->assembleL2matrices(A,B,integrand,continuous))
return false;
#if SP_DEBUG > 1
std::cout <<"---- Matrix A -----\n"<< A
<<"-------------------"<< std::endl;
std::cout <<"---- Vector B -----"<< B
<<"-------------------"<< std::endl;
#endif
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
for (i = 1; i <= nnod; i++)
for (j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
#if SP_DEBUG > 1
std::cout <<"- Solution Vector -"<< sField
<<"-------------------"<< std::endl;
#endif
return true;
}

17
src/ASM/LR/ASMu2D.h
View File

@ -377,14 +377,6 @@ public:
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool regular = true) const;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const;
//! \brief Projects inhomogenuous (scalar) dirichlet conditions by continuous L2-fit.
//! \param[in] edge low-level edge information needed to do integration
//! \param[in] values inhomogenuous function which is to be fitted
@ -433,6 +425,15 @@ protected:
// Internal utility methods
// ========================
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const;
//! \brief Connects all matching nodes on two adjacent boundary edges.
//! \param[in] edge Local edge index of this patch, in range [1,4]
//! \param neighbor The neighbor patch

29
src/ASM/LR/ASMu2Dmx.C
View File

@ -38,18 +38,18 @@
#include <numeric>
ASMu2Dmx::ASMu2Dmx (unsigned char n_s,
const std::vector<unsigned char>& n_f)
ASMu2Dmx::ASMu2Dmx (unsigned char n_s, const CharVec& n_f)
: ASMu2D(n_s), ASMmxBase(n_f)
{
threadBasis = nullptr;
}
ASMu2Dmx::ASMu2Dmx (const ASMu2Dmx& patch,
const std::vector<unsigned char>& n_f)
ASMu2Dmx::ASMu2Dmx (const ASMu2Dmx& patch, const CharVec& n_f)
: ASMu2D(patch), ASMmxBase(n_f[0]==0?patch.nfx:n_f)
{
m_basis = patch.m_basis;
threadBasis = patch.threadBasis;
nfx = patch.nfx;
nb = patch.nb;
}
@ -270,7 +270,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
if (!ok)
continue;
int iel = threadGroups[0][t][e] + 1;
auto el1 = this->getBasis(threadBasis)->elementBegin()+iel-1;
auto el1 = threadBasis->elementBegin()+iel-1;
double uh = ((*el1)->umin()+(*el1)->umax())/2.0;
double vh = ((*el1)->vmin()+(*el1)->vmax())/2.0;
std::vector<size_t> els;
@ -725,7 +725,7 @@ Vec3 ASMu2Dmx::getCoord (size_t inod) const
{
size_t b = 0;
size_t nbb = 0;
while (nbb + nb[b] < inod && b < nb.size())
while (b < nb.size() && nbb+nb[b] < inod)
nbb += nb[b++];
++b;
@ -745,20 +745,19 @@ void ASMu2Dmx::generateThreadGroups (const Integrand& integrand, bool silence,
{
// TODO: Support for div-compatible
int p1 = 0;
for (size_t i = 1; i <= m_basis.size(); ++i) {
if (this->getBasis(i)->order(0) > p1)
p1 = this->getBasis(i)->order(0), threadBasis = i;
}
for (size_t i = 1; i <= m_basis.size(); ++i)
if (this->getBasis(i)->order(0) > p1) {
threadBasis = this->getBasis(i);
p1 = threadBasis->order(0);
}
LR::generateThreadGroups(threadGroups, this->getBasis(threadBasis));
LR::generateThreadGroups(threadGroups,threadBasis);
if (silence || threadGroups[0].size() < 2) return;
std::cout <<"\nMultiple threads are utilized during element assembly.";
for (size_t i = 0; i < threadGroups[0].size(); i++)
{
std::cout <<"\n Color "<< i+1;
std::cout << ": "<< threadGroups[0][i].size() <<" elements";
}
std::cout <<"\n Color "<< i+1 <<": "
<< threadGroups[0][i].size() <<" elements";
}

15
src/ASM/LR/ASMu2Dmx.h
View File

@ -183,17 +183,26 @@ public:
int = 0, bool coordCheck = true, int thick = 1);
protected:
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const;
using ASMu2D::generateThreadGroups;
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] silence If \e true, suppress threading group outprint
//! \param[in] ignoreGlobalLM If \e true ignore global multipliers in sanity check
//! \param[in] ignoreGlobalLM If \e true, ignore global multipliers in sanity check
void generateThreadGroups(const Integrand& integrand, bool silence,
bool ignoreGlobalLM);
private:
std::vector<std::shared_ptr<LR::LRSplineSurface>> m_basis;
int threadBasis = 1; //!< Basis for thread groups
std::vector<std::shared_ptr<LR::LRSplineSurface>> m_basis; //!< All bases
LR::LRSplineSurface* threadBasis; //!< Basis for thread groups
};
#endif

75
src/ASM/LR/ASMu2Dmxrecovery.C
View File

@ -15,18 +15,19 @@
#include "LRSpline/Element.h"
#include "ASMu2Dmx.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
#include "SparseMatrix.h"
#include "DenseMatrix.h"
#include "SplineUtils.h"
#include "Utilities.h"
#include "Profiler.h"
#include "IntegrandBase.h"
#include <numeric>
/*!
\brief Expands a tensor parametrization point to an unstructured one.
*/
static void expandTensorGrid (const RealArray* in, RealArray* out)
{
out[0].resize(in[0].size()*in[1].size());
@ -41,14 +42,10 @@ static void expandTensorGrid (const RealArray* in, RealArray* out)
}
bool ASMu2Dmx::globalL2projection (Matrix& sField,
bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
if (!m_basis[0] || !m_basis[1]) return true; // silently ignore empty patches
PROFILE2("ASMu2Dmx::globalL2");
const int p1 = m_basis[0]->order(0);
const int p2 = m_basis[0]->order(1);
@ -61,17 +58,10 @@ bool ASMu2Dmx::globalL2projection (Matrix& sField,
if (!xg || !yg) return false;
if (continuous && !wg) return false;
// Set up the projection matrices
const size_t nnod = std::inner_product(nb.begin(), nb.end(), nfx.begin(), 0);
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod);
A.redim(nnod,nnod);
double dA = 0.0;
Vectors phi(m_basis.size());
std::vector<Matrix> dNdu(m_basis.size());
Matrix Xnod, Jac;
Matrices dNdu(m_basis.size());
Matrix sField, Xnod, Jac;
std::vector<Go::BasisDerivsSf> spl1(m_basis.size());
std::vector<Go::BasisPtsSf> spl0(m_basis.size());
@ -156,9 +146,8 @@ bool ASMu2Dmx::globalL2projection (Matrix& sField,
for (size_t jj = 0; jj < phi[b].size(); jj++)
{
int jnod = MNPC[iel-1][jj+el_ofs]+1;
for (size_t k=1;k<=nfx[b];++k) {
for (size_t k=1;k<=nfx[b];++k)
A((inod-nod_ofs)*nfx[b]+k+eq_ofs,(jnod-nod_ofs)*nfx[b]+k+eq_ofs) += phi[b][ii]*phi[b][jj]*dJw;
}
}
for (size_t k=1;k<=nfx[b];++k)
B((inod-nod_ofs)*nfx[b]+k+eq_ofs) += phi[b][ii]*sField(k,ip+1)*dJw;
@ -170,10 +159,32 @@ bool ASMu2Dmx::globalL2projection (Matrix& sField,
}
}
#if SP_DEBUG > 2
std::cout <<"---- Matrix A -----"<< A
return true;
}
bool ASMu2Dmx::globalL2projection (Matrix& sField,
const IntegrandBase& integrand,
bool continuous) const
{
for (size_t b = 0; b < m_basis.size(); b++)
if (!m_basis[b]) return true; // silently ignore empty patches
PROFILE2("ASMu2Dmx::globalL2");
// Assemble the projection matrices
size_t nnod = std::inner_product(nb.begin(), nb.end(), nfx.begin(), 0);
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod);
A.redim(nnod,nnod);
if (!this->assembleL2matrices(A,B,integrand,continuous))
return false;
#if SP_DEBUG > 1
std::cout <<"---- Matrix A -----\n"<< A
<<"-------------------"<< std::endl;
std::cout <<"---- Vector B -----"<< B
std::cout <<"---- Vector B -----\n"<< B
<<"-------------------"<< std::endl;
#endif
@ -184,15 +195,15 @@ bool ASMu2Dmx::globalL2projection (Matrix& sField,
sField.resize(*std::max_element(nfx.begin(), nfx.end()),
std::accumulate(nb.begin(), nb.end(), 0));
size_t eq_ofs = 0;
size_t ofs = 0;
for (size_t b = 0; b < m_basis.size(); ++b) {
for (size_t i = 1; i <= nb[b]; i++)
for (size_t j = 1; j <= nfx[b]; j++)
sField(j,i+ofs) = B((i-1)*nfx[b]+j+eq_ofs);
eq_ofs += nb[b]*nfx[b];
ofs += nb[b];
}
size_t i, j, inod = 1, jnod = 1;
for (size_t b = 0; b < m_basis.size(); b++)
for (i = 1; i <= nb[b]; i++, inod++)
for (j = 1; j <= nfx[b]; j++, jnod++)
sField(j,inod) = B(jnod);
#if SP_DEBUG > 1
std::cout <<"- Solution Vector -"<< sField
<<"-------------------"<< std::endl;
#endif
return true;
}

33
src/ASM/LR/ASMu2Drecovery.C
View File

@ -98,13 +98,11 @@ LR::LRSpline* ASMu2D::evalSolution (const IntegrandBase& integrand) const
}
bool ASMu2D::globalL2projection (Matrix& sField,
bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
if (!lrspline) return true; // silently ignore empty patches
PROFILE2("ASMu2D::globalL2");
const size_t nnod = this->getNoNodes();
const int p1 = lrspline->order(0);
const int p2 = lrspline->order(1);
@ -118,14 +116,6 @@ bool ASMu2D::globalL2projection (Matrix& sField,
if (!xg || !yg) return false;
if (continuous && !wg) return false;
// Set up the projection matrices
const size_t nnod = this->getNoNodes();
const size_t ncomp = integrand.getNoFields();
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod*ncomp);
A.redim(nnod,nnod);
double dA = 0.0;
Vector phi;
Matrix dNdu, Xnod, Jac;
@ -156,6 +146,7 @@ bool ASMu2D::globalL2projection (Matrix& sField,
expandTensorGrid(gpar.data(),unstrGpar.data());
// Evaluate the secondary solution at all integration points
Matrix sField;
if (!this->evalSolution(sField,integrand,unstrGpar.data()))
return false;
@ -195,28 +186,12 @@ bool ASMu2D::globalL2projection (Matrix& sField,
int jnod = MNPC[iel-1][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= ncomp; r++)
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
}
}
#if SP_DEBUG > 2
std::cout <<"---- Matrix A -----"<< A
<<"-------------------"<< std::endl;
std::cout <<"---- Vector B -----"<< B
<<"-------------------"<< std::endl;
#endif
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
for (size_t i = 1; i <= nnod; i++)
for (size_t j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
return true;
}

18
src/ASM/LR/ASMu3D.h
View File

@ -355,14 +355,6 @@ public:
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool regular = true) const;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const;
using ASMunstruct::generateThreadGroups;
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] integrand Object with problem-specific data and methods
@ -372,9 +364,19 @@ public:
bool ignoreGlobalLM);
protected:
// Internal utility methods
// ========================
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e false, a discrete L2-projection is used
virtual bool assembleL2matrices(SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const;
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch

33
src/ASM/LR/ASMu3Drecovery.C
View File

@ -100,13 +100,11 @@ LR::LRSpline* ASMu3D::evalSolution (const IntegrandBase& integrand) const
}
bool ASMu3D::globalL2projection (Matrix& sField,
bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
if (!lrspline) return true; // silently ignore empty patches
PROFILE2("ASMu3D::globalL2");
const size_t nnod = this->getNoNodes();
const int p1 = lrspline->order(0);
const int p2 = lrspline->order(1);
@ -123,14 +121,6 @@ bool ASMu3D::globalL2projection (Matrix& sField,
if (!xg || !yg || !zg) return false;
if (continuous && !wg) return false;
// Set up the projection matrices
const size_t nnod = this->getNoNodes();
const size_t ncomp = integrand.getNoFields();
SparseMatrix A(SparseMatrix::SUPERLU);
StdVector B(nnod*ncomp);
A.redim(nnod,nnod);
double dA = 0.0;
Vector phi;
Matrix dNdu, Xnod, Jac;
@ -162,6 +152,7 @@ bool ASMu3D::globalL2projection (Matrix& sField,
expandTensorGrid(gpar.data(),unstrGpar.data());
// Evaluate the secondary solution at all integration points
Matrix sField;
if (!this->evalSolution(sField,integrand,unstrGpar.data()))
return false;
@ -202,28 +193,12 @@ bool ASMu3D::globalL2projection (Matrix& sField,
int jnod = MNPC[iel-1][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= ncomp; r++)
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
}
}
#if SP_DEBUG > 2
std::cout <<"---- Matrix A -----"<< A
<<"-------------------"<< std::endl;
std::cout <<"---- Vector B -----"<< B
<<"-------------------"<< std::endl;
#endif
// Solve the patch-global equation system
if (!A.solve(B)) return false;
// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
for (size_t i = 1; i <= nnod; i++)
for (size_t j = 1; j <= ncomp; j++)
sField(j,i) = B(i+(j-1)*nnod);
return true;
}