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

changed: introduce a dedicated projection basis in ASMxxDmx

Browse Source 
implement ASMs2Dmx::assembleL2Matrices
implement ASMu3Dmx::assembleL2Matrices
This commit is contained in:
Arne Morten Kvarving
2017-11-06 11:25:48 +01:00
parent 85f392de13
commit d1dca5088c
16 changed files with 649 additions and 195 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/ASM/ASMbase.h
View File

@@ -573,6 +573,8 @@ public:
//! \note The implementation of this method is placed in GlbL2projector.C
bool L2projection(Matrix& fVals, FunctionBase* function, double t = 0.0);
//! \brief Returns the number of projection nodes for this patch.
virtual size_t getNoProjectionNodes() const { return this->getNoNodes(1); }
// Methods for result extraction
// =============================

17
src/ASM/ASMs2Dmx.C
View File

@@ -175,8 +175,17 @@ bool ASMs2Dmx::generateFEMTopology ()
{
if (!surf) return false;
if (m_basis.empty())
if (m_basis.empty()) {
m_basis = ASMmxBase::establishBases(surf, ASMmxBase::Type);
// we need to project on something that is not one of our bases
if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 ||
ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE)
projBasis = ASMmxBase::establishBases(surf,
ASMmxBase::FULL_CONT_RAISE_BASIS1).front();
else
projBasis = m_basis.front();
}
delete surf;
geo = surf = m_basis[geoBasis-1]->clone();
@@ -1206,3 +1215,9 @@ void ASMs2Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
for (size_t b = 1; b <= this->getNoBasis(); ++b)
this->ASMs2D::getBoundaryNodes(lIndex, nodes, b, thick, 0, local);
}
size_t ASMs2Dmx::getNoProjectionNodes() const
{
return projBasis->numCoefs_u() * projBasis->numCoefs_v();
}

4
src/ASM/ASMs2Dmx.h
View File

@@ -206,6 +206,9 @@ public:
virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
unsigned char = 0, int basis = 0) const;
//! \brief Returns the number of projection nodes for this patch.
virtual size_t getNoProjectionNodes() const;
using ASMs2D::generateThreadGroups;
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] integrand Object with problem-specific data and methods
@@ -241,6 +244,7 @@ protected:
bool continuous) const;
std::vector<std::shared_ptr<Go::SplineSurface>> m_basis; //!< Vector of bases
std::shared_ptr<Go::SplineSurface> projBasis; //!< Basis to project onto
};
#endif

154
src/ASM/ASMs2Dmxrecovery.C Normal file
View File

@@ -0,0 +1,154 @@
// $Id$
//==============================================================================
//!
//! \file ASMs2Dmxrecovery.C
//!
//! \date Dec 13 2017
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Recovery of secondary solutions for structured 2D mixed spline FE models.
//!
//==============================================================================
#include "GoTools/geometry/SplineSurface.h"
#include "ASMs2Dmx.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
#include "SparseMatrix.h"
#include "SplineUtils.h"
#include "Utilities.h"
#include <array>
bool ASMs2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
const size_t nnod = projBasis->numCoefs_u()*projBasis->numCoefs_v();
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int p11 = projBasis->order_u();
const int p21 = projBasis->order_v();
const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v();
const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1;
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr;
if (!xg || !yg) return false;
if (continuous && !wg) return false;
// Compute parameter values of the Gauss points over the whole patch
Matrix gp;
std::array<RealArray,2> gpar;
gpar[0] = this->getGaussPointParameters(gp,0,ng1,xg);
gpar[1] = this->getGaussPointParameters(gp,1,ng2,yg);
// Evaluate basis functions at all integration points
std::vector<Go::BasisPtsSf> spl0;
std::array<std::vector<Go::BasisDerivsSf>,2> spl1;
if (continuous) {
projBasis->computeBasisGrid(gpar[0],gpar[1],spl1[0]);
surf->computeBasisGrid(gpar[0],gpar[1],spl1[1]);
} else
projBasis->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 <<" *** ASMs2Dmx::assembleL2matrices: Failed for patch "<< idx+1
<<" nPoints="<< gpar[0].size()*gpar[1].size() << std::endl;
return false;
}
double dA = 1.0;
std::array<Vector, 2> phi;
phi[0].resize(p11*p21);
phi[1].resize(p1*p2);
std::array<Matrix,2> dNdu;
Matrix Xnod, J;
// === Assembly loop over all elements in the patch ==========================
int iel = 0;
for (int i2 = 0; i2 < nel2; i2++)
for (int i1 = 0; i1 < nel1; i1++, iel++)
{
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 ((dA = 0.25*this->getParametricArea(1+iel)) < 0.0)
return false; // topology error (probably logic error)
}
int ip = (i2*ng1*nel1 + i1)*ng2;
IntVec lmnpc;
if (projBasis != m_basis[0]) {
lmnpc.reserve(phi[0].size());
int vidx = (spl1[0][ip].left_idx[1]-p21+1)*projBasis->numCoefs_u();
for (int j = 0; j < p21; ++j, vidx += projBasis->numCoefs_u())
for (int i = 0; i < p11; ++i)
if (continuous)
lmnpc.push_back(spl1[0][ip].left_idx[0]-p11+1+i+vidx);
else
lmnpc.push_back(spl0[ip].left_idx[0]-p11+1+i+vidx);
}
const IntVec& mnpc = projBasis == m_basis[0] ? MNPC[iel] : lmnpc;
// --- Integration loop over all Gauss points in each direction ----------
Matrix eA(p11*p21, p11*p21);
Vectors eB(sField.rows(), Vector(p11*p21));
for (int j = 0; j < ng2; j++, ip += ng1*(nel1-1))
for (int i = 0; i < ng1; i++, ip++)
{
if (continuous) {
SplineUtils::extractBasis(spl1[0][ip],phi[0],dNdu[0]);
SplineUtils::extractBasis(spl1[1][ip],phi[1],dNdu[1]);
}
else
phi[0] = spl0[ip].basisValues;
// Compute the Jacobian inverse and derivatives
double dJw = 1.0;
if (continuous)
{
dJw = dA*wg[i]*wg[j]*utl::Jacobian(J,dNdu[1],Xnod,dNdu[1],false);
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the mass matrix
eA.outer_product(phi[0], phi[0], true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi[0],sField(r,ip+1)*dJw);
}
for (int i = 0; i < p11*p21; ++i) {
for (int j = 0; j < p11*p21; ++j)
A(mnpc[i]+1, mnpc[j]+1) += eA(i+1, j+1);
int jp = mnpc[i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;
}

107
src/ASM/ASMs2Drecovery.C
View File

@@ -605,110 +605,3 @@ Go::SplineSurface* ASMs2D::projectSolutionLocalApprox (const IntegrandBase& inte
surf->rational(),
weights);
}
bool ASMs2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
const size_t nnod = this->getNoNodes(1);
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v();
const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1;
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr;
if (!xg || !yg) return false;
if (continuous && !wg) return false;
// Compute parameter values of the Gauss points over the whole patch
Matrix gp;
std::array<RealArray,2> gpar;
gpar[0] = this->getGaussPointParameters(gp,0,ng1,xg);
gpar[1] = this->getGaussPointParameters(gp,1,ng2,yg);
// Evaluate basis functions at all integration points
std::vector<Go::BasisPtsSf> spl0, spl01;
std::vector<Go::BasisDerivsSf> spl1;
if (continuous)
surf->computeBasisGrid(gpar[0],gpar[1],spl1);
else
surf->computeBasisGrid(gpar[0],gpar[1],spl0);
m_basis[0]->computeBasisGrid(gpar[0],gpar[1],spl01);
// Evaluate the secondary solution at all integration points
Matrix sField;
if (!this->evalSolution(sField,integrand,gpar.data()))
{
std::cerr <<" *** ASMs2D::assembleL2matrices: Failed for patch "<< idx+1
<<" nPoints="<< gpar[0].size()*gpar[1].size() << std::endl;
return false;
}
double dA = 1.0;
Vector phi;
Matrix dNdu, Xnod, J;
// === Assembly loop over all elements in the patch ==========================
int iel = 0;
for (int i2 = 0; i2 < nel2; i2++)
for (int i1 = 0; i1 < nel1; i1++, iel++)
{
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 ((dA = 0.25*this->getParametricArea(1+iel)) < 0.0)
return false; // topology error (probably logic error)
}
// --- Integration loop over all Gauss points in each direction ----------
int ip = (i2*ng1*nel1 + i1)*ng2;
for (int j = 0; j < ng2; j++, ip += ng1*(nel1-1))
for (int i = 0; i < ng1; i++, ip++)
{
if (continuous)
SplineUtils::extractBasis(spl1[ip],phi,dNdu);
phi = spl01[ip].basisValues;
// Compute the Jacobian inverse and derivatives
double dJw = 1.0;
if (continuous)
{
dJw = dA*wg[i]*wg[j]*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++)
{
int inod = MNPC[iel][ii]+1;
for (size_t jj = 0; jj < phi.size(); jj++)
{
int jnod = MNPC[iel][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
}
}
return true;
}

19
src/ASM/ASMs3Dmx.C
View File

@@ -175,9 +175,18 @@ bool ASMs3Dmx::generateFEMTopology ()
{
if (!svol) return false;
if (m_basis.empty())
if (m_basis.empty()) {
m_basis = ASMmxBase::establishBases(svol, ASMmxBase::Type);
// we need to project on something that is not one of our bases
if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 ||
ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE)
projBasis = ASMmxBase::establishBases(svol,
ASMmxBase::FULL_CONT_RAISE_BASIS1).front();
else
projBasis = m_basis.front();
}
delete svol;
geo = svol = m_basis[geoBasis-1]->clone();
@@ -1348,3 +1357,11 @@ void ASMs3Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
for (size_t b = 1; b <= this->getNoBasis(); ++b)
this->ASMs3D::getBoundaryNodes(lIndex, nodes, b, thick, 0, local);
}
size_t ASMs3Dmx::getNoProjectionNodes() const
{
return projBasis->numCoefs(0) *
projBasis->numCoefs(1) *
projBasis->numCoefs(2);
}

14
src/ASM/ASMs3Dmx.h
View File

@@ -198,6 +198,9 @@ public:
virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
unsigned char = 0, int basis = 0) const;
//! \brief Returns the number of projection nodes for this patch.
virtual size_t getNoProjectionNodes() const;
//! \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
@@ -233,8 +236,17 @@ protected:
virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis = 0,
int thick = 1, int = 0, bool local = false) const;
private:
//! \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;
std::vector<std::shared_ptr<Go::SplineVolume>> m_basis; //!< Vector of bases
std::shared_ptr<Go::SplineVolume> projBasis; //!< Basis to project onto
};
#endif

170
src/ASM/ASMs3Dmxrecovery.C Normal file
View File

@@ -0,0 +1,170 @@
// $Id$
//==============================================================================
//!
//! \file ASMs3Dmxrecovery.C
//!
//! \date Dec 13 2017
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Recovery of secondary solutions for structured 3D mixed spline FE models.
//!
//==============================================================================
#include "GoTools/trivariate/SplineVolume.h"
#include "ASMs3Dmx.h"
#include "FiniteElement.h"
#include "Field.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
#include "SparseMatrix.h"
#include "SplineUtils.h"
#include "Utilities.h"
#include <array>
bool ASMs3Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
const size_t nnod = projBasis->numCoefs(0) *
projBasis->numCoefs(1) *
projBasis->numCoefs(2);
const int p1 = svol->order(0);
const int p2 = svol->order(1);
const int p3 = svol->order(2);
const int p11 = projBasis->order(0);
const int p21 = projBasis->order(1);
const int p31 = projBasis->order(2);
const int n1 = svol->numCoefs(0);
const int n2 = svol->numCoefs(1);
const int n3 = svol->numCoefs(2);
const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1;
const int nel3 = n3 - p3 + 1;
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1;
const int ng3 = continuous ? nGauss : p3 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2);
const double* zg = GaussQuadrature::getCoord(ng3);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : 0;
if (!xg || !yg || !zg) return false;
if (continuous && !wg) return false;
// Compute parameter values of the Gauss points over the whole patch
Matrix gp;
std::array<RealArray,3> gpar;
gpar[0] = this->getGaussPointParameters(gp,0,ng1,xg);
gpar[1] = this->getGaussPointParameters(gp,1,ng2,yg);
gpar[2] = this->getGaussPointParameters(gp,2,ng3,zg);
// Evaluate basis functions at all integration points
std::vector<Go::BasisPts> spl0;
std::array<std::vector<Go::BasisDerivs>,2> spl1;
if (continuous) {
projBasis->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl1[0]);
svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl1[1]);
} else
projBasis->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::assembleL2matrices: Failed for patch "<< idx+1
<<" nPoints="<< gpar[0].size()*gpar[1].size()*gpar[2].size()
<< std::endl;
return false;
}
double dV = 1.0;
std::array<Vector,2> phi;
phi[0].resize(p11*p21*p31);
phi[1].resize(p1*p2*p3);
std::array<Matrix,2> dNdu;
Matrix Xnod, J;
// === Assembly loop over all elements in the patch ==========================
int iel = 0;
for (int i3 = 0; i3 < nel3; i3++)
for (int i2 = 0; i2 < nel2; i2++)
for (int i1 = 0; i1 < nel1; i1++, iel++)
{
if (MLGE[iel] < 1) continue; // zero-volume element
if (continuous)
{
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,1+iel))
return false;
else if ((dV = 0.125*this->getParametricVolume(1+iel)) < 0.0)
return false; // topology error (probably logic error)
}
int ip = ((i3*ng2*nel2 + i2)*ng1*nel1 + i1)*ng3;
IntVec lmnpc;
if (projBasis != m_basis[0]) {
lmnpc.reserve(phi[0].size());
int nuv = projBasis->numCoefs(0)*projBasis->numCoefs(1);
int widx = (spl1[0][ip].left_idx[2]-p31+1)*nuv;
for (int k = 0; k < p31; ++k, widx += nuv) {
int vidx = (spl1[0][ip].left_idx[1]-p21+1)*projBasis->numCoefs(0);
for (int j = 0; j < p21; ++j, vidx += projBasis->numCoefs(0))
for (int i = 0; i < p11; ++i)
if (continuous)
lmnpc.push_back(spl1[0][ip].left_idx[0]-p11+1+i+vidx+widx);
else
lmnpc.push_back(spl0[ip].left_idx[0]-p11+1+i+vidx+widx);
}
}
const IntVec& mnpc = projBasis == m_basis[0] ? MNPC[iel] : lmnpc;
// --- Integration loop over all Gauss points in each direction --------
Matrix eA(p11*p21*p31, p11*p21*p31);
Vectors eB(sField.rows(), Vector(p11*p21*p31));
for (int k = 0; k < ng3; k++, ip += ng2*(nel2-1)*ng1*nel1)
for (int j = 0; j < ng2; j++, ip += ng1*(nel1-1))
for (int i = 0; i < ng1; i++, ip++)
{
if (continuous) {
SplineUtils::extractBasis(spl1[0][ip],phi[0],dNdu[0]);
SplineUtils::extractBasis(spl1[1][ip],phi[1],dNdu[1]);
} else
phi[0] = spl0[ip].basisValues;
// Compute the Jacobian inverse and derivatives
double dJw = dV;
if (continuous)
{
dJw *= wg[i]*wg[j]*wg[k]*utl::Jacobian(J,dNdu[1],Xnod,dNdu[1],false);
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the mass matrix
eA.outer_product(phi[0], phi[0], true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi[0],sField(r,ip+1)*dJw);
}
for (int i = 0; i < p11*p21*p31; ++i) {
for (int j = 0; j < p11*p21*p31; ++j)
A(mnpc[i]+1, mnpc[j]+1) += eA(i+1, j+1);
int jp = mnpc[i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;
}

2
src/ASM/GlbL2projector.C
View File

@@ -266,7 +266,7 @@ bool ASMbase::globalL2projection (Matrix& sField,
PROFILE2("ASMbase::globalL2");
// Assemble the projection matrices
size_t i, nnod = this->getNoNodes(1);
size_t i, nnod = this->getNoProjectionNodes();
size_t j, ncomp = integrand.getNoFields(2);
SparseMatrix* A;
StdVector* B;

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

@@ -182,6 +182,17 @@ bool ASMu2Dmx::generateFEMTopology ()
m_basis.resize(vec.size());
for (size_t i=0;i<vec.size();++i)
m_basis[i].reset(new LR::LRSplineSurface(vec[i].get()));
// we need to project on something that is not one of our bases
if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 ||
ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE ||
ASMmxBase::Type == ASMmxBase::SUBGRID) {
auto vec2 = ASMmxBase::establishBases(tensorspline,
ASMmxBase::FULL_CONT_RAISE_BASIS1);
projBasis.reset(new LR::LRSplineSurface(vec2.front().get()));
projBasis->generateIDs();
} else
projBasis = m_basis[0];
}
lrspline = m_basis[geoBasis-1];
@@ -839,7 +850,7 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool) const
{
#ifdef SP_DEBUG
std::cout <<"ASMu2D::evalSolution(Matrix&,const IntegrandBase&,const RealArray*,bool)\n";
std::cout <<"ASMu2Dmx::evalSolution(Matrix&,const IntegrandBase&,const RealArray*,bool)\n";
#endif
sField.resize(0,0);
@@ -1139,3 +1150,9 @@ void ASMu2Dmx::remapErrors(RealArray& errors,
errors[b->getId()] += origErr[gEl-1];
}
}
size_t ASMu2Dmx::getNoProjectionNodes() const
{
return projBasis->nBasisFunctions();
}

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

@@ -163,13 +163,8 @@ public:
virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
unsigned char = 0, int basis = 0) 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 Returns the number of projection nodes for this patch.
virtual size_t getNoProjectionNodes() const;
using ASMu2D::refine;
//! \brief Refines the mesh adaptively.
@@ -229,6 +224,7 @@ protected:
private:
std::vector<std::shared_ptr<LR::LRSplineSurface>> m_basis; //!< All bases
LR::LRSplineSurface* threadBasis; //!< Basis for thread groups
std::shared_ptr<LR::LRSplineSurface> projBasis; //!< Basis to project onto
};
#endif

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

@@ -46,8 +46,8 @@ bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
const int p1 = m_basis[0]->order(0);
const int p2 = m_basis[0]->order(1);
const int p1 = projBasis->order(0);
const int p2 = projBasis->order(1);
// Get Gaussian quadrature points
const int ng1 = continuous ? nGauss : p1 - 1;
@@ -58,6 +58,7 @@ bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (!xg || !yg) return false;
if (continuous && !wg) return false;
size_t nnod = this->getNoProjectionNodes();
double dA = 0.0;
std::array<Vector, 2> phi;
std::array<Matrix, 2> dNdu;
@@ -68,30 +69,27 @@ bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
// === Assembly loop over all elements in the patch ==========================
std::vector<LR::Element*>::iterator el1 = m_basis[geoBasis-1]->elementBegin();
for (int iel = 1; el1 != m_basis[geoBasis-1]->elementEnd(); ++el1, ++iel)
for (const LR::Element* el1 : m_basis[geoBasis-1]->getAllElements())
{
double uh = ((*el1)->umin()+(*el1)->umax())/2.0;
double vh = ((*el1)->vmin()+(*el1)->vmax())/2.0;
double uh = (el1->umin()+el1->umax())/2.0;
double vh = (el1->vmin()+el1->vmax())/2.0;
std::array<size_t, 2> els;
els[0] = m_basis[0]->getElementContaining(uh,vh)+1;
els[0] = projBasis->getElementContaining(uh,vh)+1;
els[1] = m_basis[geoBasis-1]->getElementContaining(uh,vh)+1;
int geoEl = els[1];
if (continuous)
{
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,geoEl))
if (!this->getElementCoordinates(Xnod,els[1]))
return false;
else if ((dA = 0.25*this->getParametricArea(geoEl)) < 0.0)
else if ((dA = 0.25*this->getParametricArea(els[1])) < 0.0)
return false; // topology error (probably logic error)
}
// Compute parameter values of the Gauss points over this element
RealArray gpar[2], unstrGpar[2];
this->getGaussPointParameters(gpar[0],0,ng1,geoEl,xg);
this->getGaussPointParameters(gpar[1],1,ng2,geoEl,yg);
this->getGaussPointParameters(gpar[0],0,ng1,els[1],xg);
this->getGaussPointParameters(gpar[1],1,ng2,els[1],yg);
// convert to unstructred mesh representation
expandTensorGrid(gpar, unstrGpar);
@@ -101,17 +99,27 @@ bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
return false;
// set up basis function size (for extractBasis subroutine)
phi[0].resize(m_basis[0]->getElement(els[0]-1)->nBasisFunctions());
phi[1].resize(m_basis[geoBasis-1]->getElement(els[1]-1)->nBasisFunctions());
const LR::Element* elm = projBasis->getElement(els[0]-1);
phi[0].resize(elm->nBasisFunctions());
phi[1].resize(el1->nBasisFunctions());
IntVec lmnpc;
if (projBasis != m_basis[0]) {
lmnpc.reserve(phi[0].size());
for (const LR::Basisfunction* f : elm->support())
lmnpc.push_back(f->getId());
}
const IntVec& mnpc = projBasis == m_basis[0] ? MNPC[els[1]-1] : lmnpc;
// --- Integration loop over all Gauss points in each direction ----------
Matrix eA(phi[0].size(), phi[0].size());
Vectors eB(sField.rows(), Vector(phi[0].size()));
int ip = 0;
for (int j = 0; j < ng2; j++)
for (int i = 0; i < ng1; i++, ip++)
{
if (continuous)
{
m_basis[0]->computeBasis(gpar[0][i], gpar[1][j], spl1[0], els[0]-1);
projBasis->computeBasis(gpar[0][i], gpar[1][j], spl1[0], els[0]-1);
SplineUtils::extractBasis(spl1[0],phi[0],dNdu[0]);
m_basis[geoBasis-1]->computeBasis(gpar[0][i], gpar[1][j],
spl1[1], els[1]-1);
@@ -119,7 +127,7 @@ bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
}
else
{
m_basis[0]->computeBasis(gpar[0][i], gpar[1][j], spl0[0], els[0]-1);
projBasis->computeBasis(gpar[0][i], gpar[1][j], spl0[0], els[0]-1);
phi[0] = spl0[0].basisValues;
}
@@ -131,66 +139,23 @@ bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the linear system A*x=B
size_t ncmp = sField.rows();
for (size_t ii = 0; ii < phi[0].size(); ii++)
{
int inod = MNPC[els[1]-1][ii];
for (size_t jj = 0; jj < phi[0].size(); jj++)
{
int jnod = MNPC[els[1]-1][jj];
for (size_t k = 1; k <= ncmp; ++k)
A(inod*ncmp+k, jnod*ncmp+k) += phi[0][ii]*phi[0][jj]*dJw;
}
for (size_t k = 1; k <= ncmp; ++k)
B(inod*ncmp+k) += phi[0][ii]*sField(k,ip+1)*dJw;
}
// Integrate the mass matrix
eA.outer_product(phi[0], phi[0], true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi[0],sField(r,ip+1)*dJw);
}
for (size_t i = 0; i < eA.cols(); ++i) {
for (size_t j = 0; j < eA.cols(); ++j)
A(mnpc[i]+1, mnpc[j]+1) += eA(i+1,j+1);
int jp = mnpc[i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
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 = integrand.getNoFields(2)*nb[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 -----\n"<< 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(integrand.getNoFields(2), nb[0]);
size_t inod = 1, jnod = 1;
for (size_t i = 1; i <= nb[0]; i++, inod++)
for (size_t j = 1; j <= integrand.getNoFields(2); j++, jnod++)
sField(j,inod) = B(jnod);
#if SP_DEBUG > 1
std::cout <<"- Solution Vector -"<< sField
<<"-------------------"<< std::endl;
#endif
return true;
}

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

@@ -179,7 +179,7 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
// Integrate the mass matrix
eA.outer_product(phi, phi, true, dJw);
// Integrate the rhs vector B
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi,sField(r,ip+1)*dJw);
}

17
src/ASM/LR/ASMu3Dmx.C
View File

@@ -178,8 +178,19 @@ bool ASMu3Dmx::generateFEMTopology ()
m_basis.resize(vec.size());
for (size_t i=0;i<vec.size();++i)
m_basis[i].reset(new LR::LRSplineVolume(vec[i].get()));
// we need to project on something that is not one of our bases
if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 ||
ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE ||
ASMmxBase::Type == ASMmxBase::SUBGRID) {
std::shared_ptr<Go::SplineVolume> otherBasis =
ASMmxBase::establishBases(tensorspline, ASMmxBase::FULL_CONT_RAISE_BASIS1).front();
projBasis.reset(new LR::LRSplineVolume(otherBasis.get()));
} else
projBasis = m_basis[0];
}
lrspline = m_basis[geoBasis-1];
projBasis->generateIDs();
nb.resize(m_basis.size());
for (size_t i=0; i < m_basis.size(); ++i)
@@ -916,3 +927,9 @@ void ASMu3Dmx::remapErrors (RealArray& errors,
errors[b->getId()] += origErr[gEl-1];
}
}
size_t ASMu3Dmx::getNoProjectionNodes() const
{
return projBasis->nBasisFunctions();
}

14
src/ASM/LR/ASMu3Dmx.h
View File

@@ -155,6 +155,9 @@ public:
virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
unsigned char = 0, int basis = 0) const;
//! \brief Returns the number of projection nodes for this patch.
virtual size_t getNoProjectionNodes() const;
using ASMu3D::refine;
//! \brief Refines the mesh adaptively.
//! \param[in] prm Input data used to control the refinement
@@ -170,8 +173,19 @@ public:
virtual void remapErrors(RealArray& errors,
const RealArray& origErr, bool elemErrors) const;
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;
private:
std::vector<std::shared_ptr<LR::LRSplineVolume>> m_basis; //!< Spline bases
std::shared_ptr<LR::LRSplineVolume> projBasis; //!< Basis to project onto
const std::vector<Matrices>& bezierExtract; //!< Bezier extraction matrices
std::vector<Matrices> myBezierExtract; //!< Bezier extraction matrices
};

178
src/ASM/LR/ASMu3Dmxrecovery.C Normal file
View File

@@ -0,0 +1,178 @@
// $Id$
//==============================================================================
//!
//! \file ASMu3Dmxrecovery.C
//!
//! \date Nov 9 2017
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Recovery techniques for unstructured mixed LR B-splines.
//!
//==============================================================================
#include "LRSpline/LRSplineVolume.h"
#include "LRSpline/Element.h"
#include "ASMu3Dmx.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
#include "SparseMatrix.h"
#include "SplineUtils.h"
#include "Profiler.h"
#include "matrix.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()*in[2].size());
out[1].resize(in[0].size()*in[1].size()*in[2].size());
out[2].resize(in[0].size()*in[1].size()*in[2].size());
size_t i, j, k, ip = 0;
for (k = 0; k < in[2].size(); k++)
for (j = 0; j < in[1].size(); j++)
for (i = 0; i < in[0].size(); i++, ip++) {
out[0][ip] = in[0][i];
out[1][ip] = in[1][j];
out[2][ip] = in[2][k];
}
}
bool ASMu3Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
const IntegrandBase& integrand,
bool continuous) const
{
const int p1 = projBasis->order(0);
const int p2 = projBasis->order(1);
const int p3 = projBasis->order(2);
// Get Gaussian quadrature points
const int ng1 = continuous ? nGauss : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1;
const int ng3 = continuous ? nGauss : p3 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2);
const double* zg = GaussQuadrature::getCoord(ng3);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr;
if (!xg || !yg || !zg) return false;
if (continuous && !wg) return false;
size_t nnod = this->getNoProjectionNodes();
double dV = 0.0;
Vectors phi(2);
Matrices dNdu(2);
Matrix sField, Xnod, Jac;
std::vector<Go::BasisDerivs> spl1(2);
std::vector<Go::BasisPts> spl0(2);
// === Assembly loop over all elements in the patch ==========================
LR::LRSplineVolume* geoVol;
if (m_basis[geoBasis-1]->nBasisFunctions() == projBasis->nBasisFunctions())
geoVol = m_basis[geoBasis-1].get();
else
geoVol = projBasis.get();
for (const LR::Element* el1 : geoVol->getAllElements())
{
double uh = (el1->umin()+el1->umax())/2.0;
double vh = (el1->vmin()+el1->vmax())/2.0;
double wh = (el1->wmin()+el1->wmax())/2.0;
std::vector<size_t> els;
els.push_back(projBasis->getElementContaining(uh, vh, wh) + 1);
els.push_back(m_basis[geoBasis-1]->getElementContaining(uh, vh, wh) + 1);
if (continuous)
{
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,els[1]))
return false;
else if ((dV = 0.25*this->getParametricVolume(els[1])) < 0.0)
return false; // topology error (probably logic error)
}
// Compute parameter values of the Gauss points over this element
RealArray gpar[3], unstrGpar[3];
this->getGaussPointParameters(gpar[0],0,ng1,els[1],xg);
this->getGaussPointParameters(gpar[1],1,ng2,els[1],yg);
this->getGaussPointParameters(gpar[2],2,ng3,els[1],zg);
// convert to unstructred mesh representation
expandTensorGrid(gpar, unstrGpar);
// Evaluate the secondary solution at all integration points
if (!this->evalSolution(sField,integrand,unstrGpar))
return false;
// set up basis function size (for extractBasis subroutine)
const LR::Element* elm = projBasis->getElement(els[0]-1);
phi[0].resize(elm->nBasisFunctions());
phi[1].resize(el1->nBasisFunctions());
IntVec lmnpc;
if (projBasis != m_basis[0]) {
lmnpc.reserve(phi[0].size());
for (const LR::Basisfunction* f : elm->support())
lmnpc.push_back(f->getId());
}
const IntVec& mnpc = projBasis == m_basis[0] ? MNPC[els[1]-1] : lmnpc;
// --- Integration loop over all Gauss points in each direction ----------
Matrix eA(phi[0].size(), phi[0].size());
Vectors eB(sField.rows(), Vector(phi[0].size()));
int ip = 0;
for (int k = 0; k < ng3; k++)
for (int j = 0; j < ng2; j++)
for (int i = 0; i < ng1; i++, ip++)
{
if (continuous)
{
projBasis->computeBasis(gpar[0][i], gpar[1][j], gpar[2][k],
spl1[0], els[0]-1);
SplineUtils::extractBasis(spl1[0],phi[0],dNdu[0]);
m_basis[geoBasis-1]->computeBasis(gpar[0][i], gpar[1][j], gpar[2][k],
spl1[1], els[1]-1);
SplineUtils::extractBasis(spl1[1], phi[1], dNdu[1]);
}
else
{
projBasis->computeBasis(gpar[0][i], gpar[1][j], gpar[2][k],
spl0[0], els[0]-1);
phi[0] = spl0[0].basisValues;
}
// Compute the Jacobian inverse and derivatives
double dJw = 1.0;
if (continuous)
{
dJw = dV*wg[i]*wg[j]*wg[k]*utl::Jacobian(Jac,dNdu[1],Xnod,dNdu[1],false);
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the mass matrix
eA.outer_product(phi[0], phi[0], true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi[0],sField(r,ip+1)*dJw);
}
for (size_t i = 0; i < eA.rows(); ++i) {
for (size_t j = 0; j < eA.cols(); ++j)
A(mnpc[i]+1, mnpc[j]+1) += eA(i+1,j+1);
int jp = mnpc[i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;
}