Unify ASMs2D(mx)::assembleL2matrices and ASMu2D(mx)::assembleL2matrices
This commit is contained in:
Arne Morten Kvarving
Knut Morten Okstad
parent
99796c2279
commit
1f3be6906f
@ -532,7 +532,7 @@ protected:
|
||||
|
||||
//! \brief Returns the area in the parameter space for an element.
|
||||
//! \param[in] iel Element index
|
||||
double getParametricArea(int iel) const;
|
||||
virtual double getParametricArea(int iel) const;
|
||||
//! \brief Returns boundary edge length in the parameter space for an element.
|
||||
//! \param[in] iel Element index
|
||||
//! \param[in] dir Local index of the boundary edge
|
||||
|
||||
@ -87,7 +87,7 @@ public:
|
||||
virtual char getNodeType(size_t inod) const;
|
||||
//! \brief Returns the area in the parameter space for an element.
|
||||
//! \param[in] iel Element index
|
||||
double getParametricArea(int iel) const;
|
||||
virtual double getParametricArea(int iel) const;
|
||||
//! \brief Returns boundary edge length in the parameter space for an element.
|
||||
//! \param[in] iel Element index
|
||||
//! \param[in] dir Local index of the boundary edge
|
||||
@ -239,15 +239,6 @@ public:
|
||||
int thick = 1, int = 0, bool local = false) 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;
|
||||
|
||||
std::vector<std::shared_ptr<Go::SplineSurface>> m_basis; //!< Vector of bases
|
||||
};
|
||||
|
||||
|
||||
@ -1,154 +0,0 @@
|
||||
// $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 = this->getNoProjectionNodes();
|
||||
|
||||
const int p1 = surf->order_u();
|
||||
const int p2 = surf->order_v();
|
||||
const int p11 = proj->order_u();
|
||||
const int p21 = proj->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) {
|
||||
proj->computeBasisGrid(gpar[0],gpar[1],spl1[0]);
|
||||
surf->computeBasisGrid(gpar[0],gpar[1],spl1[1]);
|
||||
} else
|
||||
proj->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 (proj != m_basis.front().get()) {
|
||||
lmnpc.reserve(phi[0].size());
|
||||
int vidx = (spl1[0][ip].left_idx[1]-p21+1)*proj->numCoefs_u();
|
||||
for (int j = 0; j < p21; ++j, vidx += proj->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 = proj == m_basis.front().get() ? 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;
|
||||
}
|
||||
@ -164,14 +164,15 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
const IntegrandBase& integrand,
|
||||
bool continuous) const
|
||||
{
|
||||
const size_t nnod = this->getNoNodes(1);
|
||||
const size_t nnod = this->getNoProjectionNodes();
|
||||
|
||||
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;
|
||||
const int g1 = surf->order_u();
|
||||
const int g2 = surf->order_v();
|
||||
const int p1 = proj->order_u();
|
||||
const int p2 = proj->order_v();
|
||||
const int n1 = proj->numCoefs_u();
|
||||
const int nel1 = surf->numCoefs_u() - g1 + 1;
|
||||
const int nel2 = surf->numCoefs_v() - g2 + 1;
|
||||
const int pmax = p1 > p2 ? p1 : p2;
|
||||
|
||||
// Get Gaussian quadrature point coordinates (and weights if continuous)
|
||||
@ -191,23 +192,26 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
|
||||
// Evaluate basis functions at all integration points
|
||||
std::vector<Go::BasisPtsSf> spl0;
|
||||
std::vector<Go::BasisDerivsSf> spl1;
|
||||
std::vector<Go::BasisDerivsSf> spl1, spl2;
|
||||
if (continuous)
|
||||
surf->computeBasisGrid(gpar[0],gpar[1],spl1);
|
||||
{
|
||||
proj->computeBasisGrid(gpar[0],gpar[1],spl1);
|
||||
surf->computeBasisGrid(gpar[0],gpar[1],spl2);
|
||||
}
|
||||
else
|
||||
surf->computeBasisGrid(gpar[0],gpar[1],spl0);
|
||||
proj->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::assembleL2matrices: Failed for patch "<< idx+1
|
||||
<<" nPoints="<< gpar[0].size()*gpar[1].size() << std::endl;
|
||||
<<" nPoints="<< gpar[0].size()*gpar[1].size() << std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
double dA = 1.0;
|
||||
Vector phi(p1*p2);
|
||||
Vector phi(p1*p2), phi2(g1*g2);
|
||||
Matrix dNdu, Xnod, J;
|
||||
|
||||
|
||||
@ -221,47 +225,66 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
|
||||
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)
|
||||
// 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 (proj != surf)
|
||||
{
|
||||
// Establish nodal point correspondance for the projection element
|
||||
int i, j, vidx;
|
||||
lmnpc.reserve(phi.size());
|
||||
if (continuous)
|
||||
vidx = (spl1[ip].left_idx[1]-p1+1)*n1 + (spl1[ip].left_idx[0]-p1+1);
|
||||
else
|
||||
vidx = (spl0[ip].left_idx[1]-p1+1)*n1 + (spl0[ip].left_idx[0]-p1+1);
|
||||
for (j = 0; j < p2; j++, vidx += n1)
|
||||
for (i = 0; i < p1; i++)
|
||||
lmnpc.push_back(vidx+i);
|
||||
}
|
||||
const IntVec& mnpc = proj == surf ? MNPC[iel] : lmnpc;
|
||||
|
||||
// --- Integration loop over all Gauss points in each direction ----------
|
||||
|
||||
Matrix eA(p1*p2, p1*p2);
|
||||
Vectors eB(sField.rows(), Vector(p1*p1));
|
||||
int ip = (i2*ng1*nel1 + i1)*ng2;
|
||||
Vectors eB(sField.rows(), Vector(p1*p2));
|
||||
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);
|
||||
else
|
||||
phi = spl0[ip].basisValues;
|
||||
for (int i = 0; i < ng1; i++, ip++)
|
||||
{
|
||||
if (continuous)
|
||||
{
|
||||
SplineUtils::extractBasis(spl1[ip],phi,dNdu);
|
||||
SplineUtils::extractBasis(spl2[ip],phi2,dNdu);
|
||||
}
|
||||
else
|
||||
phi = 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,Xnod,dNdu,false);
|
||||
if (dJw == 0.0) continue; // skip singular points
|
||||
}
|
||||
// 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 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);
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < p1*p2; ++i) {
|
||||
for (int j = 0; j < p1*p2; ++j)
|
||||
A(MNPC[iel][i]+1, MNPC[iel][j]+1) += eA(i+1, j+1);
|
||||
A(mnpc[i]+1, mnpc[j]+1) += eA(i+1, j+1);
|
||||
|
||||
int jp = MNPC[iel][i]+1;
|
||||
int jp = mnpc[i]+1;
|
||||
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
|
||||
B(jp) += eB[r](1+i);
|
||||
}
|
||||
|
||||
@ -215,15 +215,6 @@ public:
|
||||
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;
|
||||
|
||||
using ASMu2D::generateThreadGroups;
|
||||
//! \brief Generates element groups for multi-threading of interior integrals.
|
||||
//! \param[in] integrand Object with problem-specific data and methods
|
||||
|
||||
@ -1,161 +0,0 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file ASMu2Dmxrecovery.C
|
||||
//!
|
||||
//! \date May 11 2015
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Recovery techniques for unstructured mixed LR B-splines.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "LRSpline/LRSplineSurface.h"
|
||||
#include "LRSpline/Element.h"
|
||||
|
||||
#include "ASMu2Dmx.h"
|
||||
#include "IntegrandBase.h"
|
||||
#include "CoordinateMapping.h"
|
||||
#include "GaussQuadrature.h"
|
||||
#include "SparseMatrix.h"
|
||||
#include "SplineUtils.h"
|
||||
#include "Profiler.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());
|
||||
out[1].resize(in[0].size()*in[1].size());
|
||||
|
||||
size_t i, j, ip = 0;
|
||||
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];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bool ASMu2Dmx::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
const IntegrandBase& integrand,
|
||||
bool continuous) const
|
||||
{
|
||||
const int p1 = projBasis->order(0);
|
||||
const int p2 = projBasis->order(1);
|
||||
|
||||
// Get Gaussian quadrature points
|
||||
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;
|
||||
|
||||
size_t nnod = this->getNoProjectionNodes();
|
||||
double dA = 0.0;
|
||||
std::array<Vector, 2> phi;
|
||||
std::array<Matrix, 2> dNdu;
|
||||
Matrix sField, Xnod, Jac;
|
||||
std::array<Go::BasisDerivsSf, 2> spl1;
|
||||
std::array<Go::BasisPtsSf, 2> spl0;
|
||||
|
||||
|
||||
// === Assembly loop over all elements in the patch ==========================
|
||||
|
||||
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;
|
||||
std::array<size_t, 2> els;
|
||||
els[0] = projBasis->getElementContaining(uh,vh)+1;
|
||||
els[1] = m_basis[geoBasis-1]->getElementContaining(uh,vh)+1;
|
||||
|
||||
if (continuous)
|
||||
{
|
||||
// Set up control point (nodal) coordinates for current element
|
||||
if (!this->getElementCoordinates(Xnod,els[1]))
|
||||
return false;
|
||||
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,els[1],xg);
|
||||
this->getGaussPointParameters(gpar[1],1,ng2,els[1],yg);
|
||||
|
||||
// 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 j = 0; j < ng2; j++)
|
||||
for (int i = 0; i < ng1; i++, ip++)
|
||||
{
|
||||
if (continuous)
|
||||
{
|
||||
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);
|
||||
SplineUtils::extractBasis(spl1[1], phi[1], dNdu[1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
projBasis->computeBasis(gpar[0][i], gpar[1][j], spl0[0], els[0]-1);
|
||||
phi[0] = spl0[0].basisValues;
|
||||
}
|
||||
|
||||
// Compute the Jacobian inverse and derivatives
|
||||
double dJw = 1.0;
|
||||
if (continuous)
|
||||
{
|
||||
dJw = dA*wg[i]*wg[j]*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.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;
|
||||
}
|
||||
@ -31,10 +31,11 @@ bool ASMu2D::getGrevilleParameters (RealArray& prm, int dir, int basisNum) const
|
||||
if (!this->getBasis(basisNum) || dir < 0 || dir > 1) return false;
|
||||
|
||||
const LR::LRSpline* lrspline = this->getBasis(basisNum);
|
||||
|
||||
prm.clear();
|
||||
prm.reserve(lrspline->nBasisFunctions());
|
||||
|
||||
for (const LR::Basisfunction *b : lrspline->getAllBasisfunctions())
|
||||
for (const LR::Basisfunction* b : lrspline->getAllBasisfunctions())
|
||||
prm.push_back(b->getGrevilleParameter()[dir]);
|
||||
|
||||
return true;
|
||||
@ -101,10 +102,10 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
const IntegrandBase& integrand,
|
||||
bool continuous) const
|
||||
{
|
||||
const size_t nnod = this->getNoNodes();
|
||||
size_t nnod = this->getNoProjectionNodes();
|
||||
|
||||
const int p1 = lrspline->order(0);
|
||||
const int p2 = lrspline->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;
|
||||
@ -116,17 +117,21 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
if (continuous && !wg) return false;
|
||||
|
||||
double dA = 0.0;
|
||||
Vector phi;
|
||||
Vector phi, phi2;
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Go::BasisDerivsSf spl1;
|
||||
Go::BasisPtsSf spl0;
|
||||
Go::BasisDerivsSf spl1, spl2;
|
||||
|
||||
|
||||
// === Assembly loop over all elements in the patch ==========================
|
||||
|
||||
size_t iel = 1;
|
||||
for (const LR::Element* el : lrspline->getAllElements())
|
||||
for (const LR::Element* el1 : lrspline->getAllElements())
|
||||
{
|
||||
double uh = (el1->umin()+el1->umax())/2.0;
|
||||
double vh = (el1->vmin()+el1->vmax())/2.0;
|
||||
int ielp = projBasis->getElementContaining(uh,vh);
|
||||
int iel = lrspline->getElementContaining(uh,vh)+1;
|
||||
|
||||
if (continuous)
|
||||
{
|
||||
// Set up control point (nodal) coordinates for current element
|
||||
@ -140,8 +145,6 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
std::array<RealArray,2> gpar, unstrGpar;
|
||||
this->getGaussPointParameters(gpar[0],0,ng1,iel,xg);
|
||||
this->getGaussPointParameters(gpar[1],1,ng2,iel,yg);
|
||||
|
||||
// convert to unstructured mesh representation
|
||||
expandTensorGrid(gpar.data(),unstrGpar.data());
|
||||
|
||||
// Evaluate the secondary solution at all integration points
|
||||
@ -149,24 +152,37 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
if (!this->evalSolution(sField,integrand,unstrGpar.data()))
|
||||
return false;
|
||||
|
||||
// set up basis function size (for extractBasis subroutine)
|
||||
phi.resize(el->nBasisFunctions());
|
||||
// Set up basis function size (for extractBasis subroutine)
|
||||
const LR::Element* elm = projBasis->getElement(ielp);
|
||||
size_t nbf = elm->nBasisFunctions();
|
||||
|
||||
IntVec lmnpc;
|
||||
if (projBasis != lrspline)
|
||||
{
|
||||
lmnpc.reserve(nbf);
|
||||
for (const LR::Basisfunction* f : elm->support())
|
||||
lmnpc.push_back(f->getId());
|
||||
}
|
||||
const IntVec& mnpc = projBasis == lrspline ? MNPC[iel-1] : lmnpc;
|
||||
|
||||
// --- Integration loop over all Gauss points in each direction ----------
|
||||
Matrix eA(MNPC[iel-1].size(), MNPC[iel-1].size());
|
||||
Vectors eB(sField.rows(), Vector(MNPC[iel-1].size()));
|
||||
|
||||
Matrix eA(nbf, nbf);
|
||||
Vectors eB(sField.rows(), Vector(nbf));
|
||||
int ip = 0;
|
||||
for (int j = 0; j < ng2; j++)
|
||||
for (int i = 0; i < ng1; i++, ip++)
|
||||
{
|
||||
if (continuous)
|
||||
{
|
||||
lrspline->computeBasis(gpar[0][i],gpar[1][j],spl1,iel-1);
|
||||
projBasis->computeBasis(gpar[0][i],gpar[1][j],spl1,ielp);
|
||||
SplineUtils::extractBasis(spl1,phi,dNdu);
|
||||
lrspline->computeBasis(gpar[0][i],gpar[1][j],spl2,iel-1);
|
||||
SplineUtils::extractBasis(spl2,phi2,dNdu);
|
||||
}
|
||||
else
|
||||
{
|
||||
lrspline->computeBasis(gpar[0][i],gpar[1][j],spl0,iel-1);
|
||||
projBasis->computeBasis(gpar[0][i],gpar[1][j],spl0,ielp);
|
||||
phi = spl0.basisValues;
|
||||
}
|
||||
|
||||
@ -186,15 +202,14 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
|
||||
eB[r-1].add(phi,sField(r,ip+1)*dJw);
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < MNPC[iel-1].size(); ++i) {
|
||||
for (size_t j = 0; j < MNPC[iel-1].size(); ++j)
|
||||
A(MNPC[iel-1][i]+1, MNPC[iel-1][j]+1) += eA(i+1, j+1);
|
||||
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[iel-1][i]+1;
|
||||
int jp = mnpc[i]+1;
|
||||
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
|
||||
B(jp) += eB[r](1+i);
|
||||
}
|
||||
++iel;
|
||||
}
|
||||
|
||||
return true;
|
||||
@ -254,7 +269,7 @@ LR::LRSplineSurface* ASMu2D::scRecovery (const IntegrandBase& integrand) const
|
||||
size_t k, l, ip = 0;
|
||||
std::vector<LR::Element*>::const_iterator elStart, elEnd, el;
|
||||
std::vector<LR::Element*> supportElements;
|
||||
for (LR::Basisfunction *b : lrspline->getAllBasisfunctions())
|
||||
for (LR::Basisfunction* b : lrspline->getAllBasisfunctions())
|
||||
{
|
||||
#if SP_DEBUG > 2
|
||||
std::cout <<"Basis: "<< *b <<"\n ng1 ="<< ng1 <<"\n ng2 ="<< ng2
|
||||
@ -301,17 +316,14 @@ LR::LRSplineSurface* ASMu2D::scRecovery (const IntegrandBase& integrand) const
|
||||
for (el = elStart; el != elEnd; ++el)
|
||||
{
|
||||
int iel = (**el).getId()+1;
|
||||
#if SP_DEBUG > 2
|
||||
std::cout <<"Element "<< **el << std::endl;
|
||||
#endif
|
||||
|
||||
// evaluate all gauss points for this element
|
||||
// Compute parameter values of the Gauss points over this element
|
||||
std::array<RealArray,2> gaussPt, unstrGauss;
|
||||
this->getGaussPointParameters(gaussPt[0],0,ng1,iel,xg);
|
||||
this->getGaussPointParameters(gaussPt[1],1,ng2,iel,yg);
|
||||
|
||||
#if SP_DEBUG > 2
|
||||
std::cout << "Element " << **el << std::endl;
|
||||
#endif
|
||||
|
||||
// convert to unstructured mesh representation
|
||||
expandTensorGrid(gaussPt.data(),unstrGauss.data());
|
||||
|
||||
// Evaluate the secondary solution at all Gauss points
|
||||
|
||||
Loading…
Reference in New Issue
Block a user