Non-homogenous dirichlet BC for 2D LR-splines
This commit is contained in:
Kjetil Andre Johannessen
Arne Morten Kvarving
parent
4e8fa1ede4
commit
1a346dba8b
@ -131,6 +131,7 @@ void ASMu2D::clear (bool retainGeometry)
|
||||
|
||||
// Erase the FE data
|
||||
this->ASMbase::clear(retainGeometry);
|
||||
this->dirich.clear();
|
||||
}
|
||||
|
||||
|
||||
@ -615,69 +616,113 @@ std::vector<int> ASMu2D::getEdgeNodes (int edge, int basis) const
|
||||
|
||||
void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
|
||||
{
|
||||
// figure out function index offset (when using multiple basis)
|
||||
size_t ofs = 1;
|
||||
for (int i = 1; i < basis; i++)
|
||||
ofs += this->getNoNodes(i);
|
||||
|
||||
// figure out what edge we are at
|
||||
LR::parameterEdge edge;
|
||||
int c1, c2;
|
||||
switch (dir) {
|
||||
case -2:
|
||||
edge = LR::SOUTH;
|
||||
c1 = this->getCorner(-1, -1, basis);
|
||||
c2 = this->getCorner( 1, -1, basis);
|
||||
case -2: edge = LR::SOUTH; break;
|
||||
case -1: edge = LR::WEST; break;
|
||||
case 1: edge = LR::EAST; break;
|
||||
case 2: edge = LR::NORTH; break;
|
||||
default: return;
|
||||
}
|
||||
|
||||
// fetch the right basis to consider
|
||||
LR::LRSplineSurface* lr = this->getBasis(basis);
|
||||
|
||||
// get all elements and functions on this edge
|
||||
std::vector<LR::Basisfunction*> edgeFunctions;
|
||||
std::vector<LR::Element*> edgeElements;
|
||||
lr->getEdgeFunctions(edgeFunctions, edge);
|
||||
lr->getEdgeElements( edgeElements, edge);
|
||||
|
||||
// find the corners since these are not to be included in the L2-fitting
|
||||
// of the inhomogenuous dirichlet boundaries; corners are interpolatory.
|
||||
// Optimization note: loop over the "edge"-container to manually pick up the
|
||||
// end nodes. LRspine::getEdgeFunctions() does a global search.
|
||||
std::vector<LR::Basisfunction*> c1, c2;
|
||||
switch (edge)
|
||||
{
|
||||
case LR::SOUTH:
|
||||
lr->getEdgeFunctions(c1, LR::SOUTH_WEST);
|
||||
lr->getEdgeFunctions(c2, LR::SOUTH_EAST);
|
||||
break;
|
||||
case -1:
|
||||
edge = LR::WEST;
|
||||
c1 = this->getCorner(-1, -1, basis);
|
||||
c2 = this->getCorner(-1, 1, basis);
|
||||
case LR::WEST:
|
||||
lr->getEdgeFunctions(c1, LR::SOUTH_WEST);
|
||||
lr->getEdgeFunctions(c2, LR::NORTH_WEST);
|
||||
break;
|
||||
case 1:
|
||||
edge = LR::EAST;
|
||||
c1 = this->getCorner( 1, -1, basis);
|
||||
c2 = this->getCorner( 1, 1, basis);
|
||||
case LR::EAST:
|
||||
lr->getEdgeFunctions(c1, LR::NORTH_EAST);
|
||||
lr->getEdgeFunctions(c2, LR::SOUTH_EAST);
|
||||
break;
|
||||
case 2:
|
||||
edge = LR::NORTH;
|
||||
c1 = this->getCorner(-1, 1, basis);
|
||||
c2 = this->getCorner( 1, 1, basis);
|
||||
case LR::NORTH:
|
||||
lr->getEdgeFunctions(c1, LR::NORTH_WEST);
|
||||
lr->getEdgeFunctions(c2, LR::NORTH_EAST);
|
||||
break;
|
||||
default: return;
|
||||
}
|
||||
std::vector<int> nodes = this->getEdgeNodes(edge,basis);
|
||||
nodes.erase(std::find(nodes.begin(), nodes.end(), c1));
|
||||
nodes.erase(std::find(nodes.begin(), nodes.end(), c2));
|
||||
|
||||
int bcode = code;
|
||||
if (code > 0) {
|
||||
std::vector<LR::Basisfunction*> funcs;
|
||||
dirich.push_back(DirichletEdge(this->getBasis(basis)->edgeCurve(edge,funcs),
|
||||
dof,code));
|
||||
} else if (code < 0)
|
||||
bcode = -code;
|
||||
// build up the local element/node correspondence needed by the projection
|
||||
// call on this edge by ASMu2D::updateDirichlet()
|
||||
DirichletEdge de(edgeFunctions.size(), edgeElements.size(), dof, code, basis);
|
||||
de.corners[0] = c1[0]->getId();
|
||||
de.corners[1] = c2[0]->getId();
|
||||
de.edg = edge;
|
||||
de.lr = lr;
|
||||
int bcode = abs(code);
|
||||
|
||||
// Skip the first and last function if we are requesting an open boundary.
|
||||
if (!open)
|
||||
this->prescribe(c1,dof,bcode);
|
||||
|
||||
for (size_t i = 0; i < nodes.size(); i++)
|
||||
if (this->prescribe(nodes[i],dof,-code) == 0 && code > 0)
|
||||
dirich.back().nodes.push_back(std::make_pair(i+1,nodes[i]));
|
||||
|
||||
if (!open)
|
||||
this->prescribe(c2,dof,bcode);
|
||||
|
||||
if (code > 0)
|
||||
if (dirich.back().nodes.empty())
|
||||
dirich.pop_back(); // In the unlikely event of a 2-point boundary
|
||||
#if SP_DEBUG > 1
|
||||
else
|
||||
int j = 0;
|
||||
for (auto b : edgeFunctions )
|
||||
{
|
||||
std::cout <<"Non-corner boundary nodes:";
|
||||
for (size_t i = 0; i < dirich.back().nodes.size(); i++)
|
||||
std::cout <<" ("<< dirich.back().nodes[i].first
|
||||
<<","<< dirich.back().nodes[i].second
|
||||
<<")";
|
||||
std::cout <<"\nThese nodes will be subjected to Dirichlet projection"
|
||||
<< std::endl;
|
||||
de.MLGN[j++] = b->getId();
|
||||
// skip corners for open boundaries
|
||||
if (open && (b->getId() == de.corners[0] || b->getId() == de.corners[1]))
|
||||
continue;
|
||||
else
|
||||
{
|
||||
// corners are interpolated (positive 'code')
|
||||
if (b->getId() == de.corners[0] || b->getId() == de.corners[1])
|
||||
this->prescribe(b->getId()+ofs, dof, bcode);
|
||||
// inhomogenuous dirichlet conditions by function evaluation (negative 'code')
|
||||
else if (code > 0)
|
||||
this->prescribe(b->getId()+ofs, dof, -bcode);
|
||||
// (in)homogenuous constant dirichlet conditions
|
||||
else
|
||||
this->prescribe(b->getId()+ofs, dof, bcode);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
// build MLGE and MNPC matrix
|
||||
for (size_t i=0; i<edgeElements.size(); i++)
|
||||
{
|
||||
LR::Element* el = edgeElements[i];
|
||||
|
||||
// for mixed FEM models, let MLGE point to the *geometry* index
|
||||
if (de.lr != this->lrspline.get())
|
||||
{
|
||||
double umid = (el->umax() + el->umin())/2.0;
|
||||
double vmid = (el->vmax() + el->vmin())/2.0;
|
||||
de.MLGE[i] = lrspline->getElementContaining(umid, vmid);
|
||||
}
|
||||
else
|
||||
{
|
||||
de.MLGE[i] = el->getId();
|
||||
}
|
||||
for (auto b : el->support())
|
||||
{
|
||||
de.MNPC[i].push_back(-1);
|
||||
for (auto l2g : de.MLGN)
|
||||
// can't do map::find, since this works on first
|
||||
if (b->getId() == l2g.second)
|
||||
de.MNPC[i].back() = l2g.first;
|
||||
}
|
||||
}
|
||||
if (code > 0)
|
||||
dirich.push_back(de);
|
||||
}
|
||||
|
||||
|
||||
@ -1773,60 +1818,67 @@ bool ASMu2D::updateDirichlet (const std::map<int,RealFunc*>& func,
|
||||
const std::map<int,VecFunc*>& vfunc, double time,
|
||||
const std::map<int,int>* g2l)
|
||||
{
|
||||
std::map<int,RealFunc*>::const_iterator fit;
|
||||
std::map<int,VecFunc*>::const_iterator vfit;
|
||||
|
||||
std::map<int,RealFunc*>::const_iterator fit;
|
||||
std::map<int,VecFunc*>::const_iterator vfit;
|
||||
std::map<int,int>::const_iterator nit;
|
||||
std::vector<DirichletEdge>::const_iterator dit;
|
||||
std::vector<Ipair>::const_iterator nit;
|
||||
|
||||
for (size_t i = 0; i < dirich.size(); i++)
|
||||
{
|
||||
// Project the function onto the spline curve basis
|
||||
Go::SplineCurve* dcrv = 0;
|
||||
// figure out function index offset (when using multiple basis)
|
||||
size_t ofs = 0;
|
||||
for (int j = 1; j < dirich[i].basis; j++)
|
||||
ofs += this->getNoNodes(j);
|
||||
|
||||
RealArray edgeControlpoints;
|
||||
Real2DMat edgeControlmatrix;
|
||||
if ((fit = func.find(dirich[i].code)) != func.end())
|
||||
dcrv = SplineUtils::project(dirich[i].curve,*fit->second,time);
|
||||
edgeL2projection(dirich[i], *fit->second, edgeControlpoints, time);
|
||||
else if ((vfit = vfunc.find(dirich[i].code)) != vfunc.end())
|
||||
dcrv = SplineUtils::project(dirich[i].curve,*vfit->second,nf,time);
|
||||
edgeL2projection(dirich[i], *vfit->second, edgeControlmatrix, time);
|
||||
else
|
||||
{
|
||||
std::cerr <<" *** ASMu2D::updateDirichlet: Code "<< dirich[i].code
|
||||
<<" is not associated with any function."<< std::endl;
|
||||
return false;
|
||||
}
|
||||
if (!dcrv)
|
||||
if (edgeControlpoints.empty() && edgeControlmatrix.empty())
|
||||
{
|
||||
std::cerr <<" *** ASMu2D::updateDirichlet: Projection failure."
|
||||
<< std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
// Loop over the (interior) nodes (control points) of this boundary curve
|
||||
for (nit = dirich[i].nodes.begin(); nit != dirich[i].nodes.end(); ++nit)
|
||||
// Loop over the nodes of this boundary curve
|
||||
for (nit = dirich[i].MLGN.begin(); nit != dirich[i].MLGN.end(); nit++)
|
||||
{
|
||||
// skip corner nodes, since these are special cased (interpolatory)
|
||||
if (nit->second == dirich[i].corners[0] ||
|
||||
nit->second == dirich[i].corners[1])
|
||||
continue;
|
||||
for (int dofs = dirich[i].dof; dofs > 0; dofs /= 10)
|
||||
{
|
||||
int dof = dofs%10;
|
||||
// Find the constraint equation for current (node,dof)
|
||||
MPC pDOF(MLGN[nit->second-1],dof);
|
||||
MPC pDOF(MLGN[nit->second]+ofs,dof);
|
||||
MPCIter mit = mpcs.find(&pDOF);
|
||||
if (mit == mpcs.end()) continue; // probably a deleted constraint
|
||||
|
||||
// Find index to the control point value for this (node,dof) in dcrv
|
||||
RealArray::const_iterator cit = dcrv->coefs_begin();
|
||||
if (dcrv->dimension() > 1) // A vector field is specified
|
||||
cit += (nit->first-1)*dcrv->dimension() + (dof-1);
|
||||
else // A scalar field is specified at this dof
|
||||
cit += (nit->first-1);
|
||||
|
||||
// Now update the prescribed value in the constraint equation
|
||||
(*mit)->setSlaveCoeff(*cit);
|
||||
#if SP_DEBUG > 1
|
||||
if (edgeControlpoints.empty()) // vector conditions
|
||||
(*mit)->setSlaveCoeff(edgeControlmatrix[dof-1][nit->first]);
|
||||
else //scalar condition
|
||||
(*mit)->setSlaveCoeff(edgeControlpoints[nit->first]);
|
||||
#if SP_DEBUG > 1
|
||||
std::cout <<"Updated constraint: "<< **mit;
|
||||
#endif
|
||||
#endif
|
||||
}
|
||||
delete dcrv;
|
||||
}
|
||||
}
|
||||
|
||||
// The parent class method takes care of the corner nodes with direct
|
||||
// evaluation of the Dirichlet functions (since they are interpolatory)
|
||||
// evaluation of the Dirichlet functions; since they are interpolatory
|
||||
return this->ASMbase::updateDirichlet(func,vfunc,time,g2l);
|
||||
}
|
||||
|
||||
|
||||
@ -16,6 +16,7 @@
|
||||
|
||||
#include "ASMunstruct.h"
|
||||
#include "ASM2D.h"
|
||||
#include "LRSpline/LRSpline.h"
|
||||
#include <memory>
|
||||
|
||||
class FiniteElement;
|
||||
@ -316,6 +317,24 @@ public:
|
||||
const int* npe, char project = '\0') const;
|
||||
|
||||
private:
|
||||
//! \brief Struct representing an inhomogeneous Dirichlet boundary condition.
|
||||
struct DirichletEdge
|
||||
{
|
||||
LR::LRSplineSurface *lr; //!< Pointer to the right object (in case of multiple bases)
|
||||
LR::parameterEdge edg; //!< Which edge is this
|
||||
IntVec MLGE; //!< Local-to-Global Element numbers
|
||||
std::map<int,int> MLGN; //!< Local-to-Global Nodal numbers
|
||||
IntMat MNPC; //!< Matrix of Nodal-Point Correpondanse
|
||||
int dof; //!< Local DOF to constrain along the boundary
|
||||
int code; //!< Inhomogeneous Dirichlet condition code
|
||||
int basis; //!< Index to the basis used
|
||||
int corners[2];//!< Index of the two end-points of this line
|
||||
|
||||
//! \brief Default constructor.
|
||||
DirichletEdge(int numbBasis, int numbElements, int d = 0, int c = 0, int b = 1)
|
||||
: MLGE(numbElements), MNPC(numbElements), dof(d), code(c), basis(b) {}
|
||||
};
|
||||
|
||||
//! \brief Projects the secondary solution field onto the primary basis.
|
||||
//! \param[in] integrand Object with problem-specific data and methods
|
||||
LR::LRSplineSurface* projectSolution(const IntegrandBase& integrand) const;
|
||||
@ -361,6 +380,26 @@ public:
|
||||
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
|
||||
//! \param[out] result fitted value in terms of control-point values
|
||||
//! \param[in] time time used in dynamic problems
|
||||
bool edgeL2projection (const DirichletEdge& edge,
|
||||
const RealFunc& values,
|
||||
RealArray& result,
|
||||
double time) const;
|
||||
|
||||
//! \brief Projects inhomogenuous (vector) 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
|
||||
//! \param[out] result fitted value in terms of control-point values
|
||||
//! \param[in] time time used in dynamic problems
|
||||
bool edgeL2projection (const DirichletEdge& edge,
|
||||
const VecFunc& values,
|
||||
Real2DMat& result,
|
||||
double time) const;
|
||||
|
||||
//! \brief Transfers Gauss point variables from old basis to this patch.
|
||||
//! \param[in] oldBasis The LR-spline basis to transfer from
|
||||
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
|
||||
@ -438,17 +477,6 @@ public:
|
||||
typedef std::pair<int,int> Ipair; //!< Convenience type
|
||||
|
||||
protected:
|
||||
//! \brief Struct representing an inhomogeneous Dirichlet boundary condition.
|
||||
struct DirichletEdge
|
||||
{
|
||||
Go::SplineCurve* curve; //!< Pointer to spline curve for the boundary
|
||||
int dof; //!< Local DOF to constrain along the boundary
|
||||
int code; //!< Inhomogeneous Dirichlet condition code
|
||||
std::vector<Ipair> nodes; //!< Nodes subjected to projection on the boundary
|
||||
//! \brief Default constructor.
|
||||
DirichletEdge(Go::SplineCurve* sc = nullptr, int d = 0, int c = 0)
|
||||
: curve(sc), dof(d), code(c) {}
|
||||
};
|
||||
|
||||
std::shared_ptr<LR::LRSplineSurface> lrspline; //!< Pointer to the LR-spline surface object
|
||||
|
||||
|
||||
@ -23,8 +23,11 @@
|
||||
#include "Utilities.h"
|
||||
#include "Profiler.h"
|
||||
#include "IntegrandBase.h"
|
||||
#include "FiniteElement.h"
|
||||
#include <array>
|
||||
|
||||
#include <fstream>
|
||||
|
||||
|
||||
bool ASMu2D::getGrevilleParameters (RealArray& prm, int dir) const
|
||||
{
|
||||
@ -457,3 +460,343 @@ LR::LRSplineSurface* ASMu2D::regularInterpolation (const RealArray& upar,
|
||||
|
||||
return ans;
|
||||
}
|
||||
|
||||
|
||||
bool ASMu2D::edgeL2projection (const DirichletEdge& edge,
|
||||
const RealFunc& values,
|
||||
RealArray& result,
|
||||
double time) const
|
||||
{
|
||||
size_t n = edge.MLGN.size();
|
||||
SparseMatrix A(SparseMatrix::SUPERLU);
|
||||
StdVector B(n);
|
||||
A.resize(n,n);
|
||||
|
||||
// Get Gaussian quadrature points and weights
|
||||
const double* xg = GaussQuadrature::getCoord(nGauss);
|
||||
const double* wg = GaussQuadrature::getWeight(nGauss);
|
||||
if (!xg || !wg) return false;
|
||||
|
||||
// find the normal and tangent direction for the edge
|
||||
int edgeDir, t1, t2;
|
||||
switch (edge.edg)
|
||||
{ // id normal tangent
|
||||
case LR::WEST: edgeDir = -1; t1=1; t2=2; break;
|
||||
case LR::EAST: edgeDir = +1; t1=1; t2=2; break;
|
||||
case LR::SOUTH: edgeDir = -2; t1=2; t2=1; break;
|
||||
case LR::NORTH: edgeDir = +2; t1=2; t2=1; break;
|
||||
default: return false;
|
||||
}
|
||||
|
||||
std::array<Vector,2> gpar;
|
||||
for (int d = 0; d < 2; d++)
|
||||
if (-1-d == edgeDir)
|
||||
{
|
||||
gpar[d].resize(nGauss);
|
||||
gpar[d].fill(d == 0 ? lrspline->startparam(0) : lrspline->startparam(1));
|
||||
}
|
||||
else if (1+d == edgeDir)
|
||||
{
|
||||
gpar[d].resize(nGauss);
|
||||
gpar[d].fill(d == 0 ? lrspline->endparam(0) : lrspline->endparam(1));
|
||||
}
|
||||
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
Vec3 normal;
|
||||
int iGp = 0;
|
||||
|
||||
// === Assembly loop over all elements on the patch edge =====================
|
||||
for (size_t i=0; i<edge.MLGE.size(); i++) // for all edge elements
|
||||
{
|
||||
FiniteElement fe(edge.MNPC[i].size());
|
||||
fe.iel = edge.MLGE[i];
|
||||
|
||||
// Get element edge length in the parameter space
|
||||
double dS = this->getParametricLength(fe.iel+1,t1);
|
||||
if (dS < 0.0) return false; // topology error (probably logic error)
|
||||
|
||||
// Set up control point coordinates for current element
|
||||
if (!this->getElementCoordinates(Xnod,fe.iel+1)) return false;
|
||||
|
||||
// Initialize element quantities
|
||||
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
|
||||
|
||||
// Get integration gauss points over this element
|
||||
this->getGaussPointParameters(gpar[t2-1],t2-1,nGauss,fe.iel+1,xg);
|
||||
|
||||
// --- Integration loop over all Gauss points along the edge -------------
|
||||
for (int j = 0; j < nGauss; j++)
|
||||
{
|
||||
fe.iGP = iGp++; // Global integration point counter
|
||||
|
||||
// Local element coordinates and parameter values
|
||||
// of current integration point
|
||||
fe.xi = xg[j];
|
||||
fe.eta = xg[j];
|
||||
fe.u = gpar[0][j];
|
||||
fe.v = gpar[1][j];
|
||||
|
||||
// Evaluate basis function (geometry) derivatives at current integration points
|
||||
Go::BasisDerivsSf spline;
|
||||
lrspline->computeBasis(fe.u, fe.v, spline, fe.iel);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
SplineUtils::extractBasis(spline,fe.N,dNdu);
|
||||
|
||||
// Compute basis function derivatives and the edge normal
|
||||
fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2);
|
||||
if (fe.detJxW == 0.0) continue; // skip singular points
|
||||
|
||||
if (edgeDir < 0) normal *= -1.0;
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.t = time;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
fe.detJxW *= 0.5*dS*wg[j];
|
||||
|
||||
// For mixed basis, we need to compute functions separate from geometry
|
||||
if (edge.lr != lrspline.get())
|
||||
{
|
||||
// different lrspline instances enumerate elements differently
|
||||
int basis_el = edge.lr->getElementContaining(fe.u, fe.v);
|
||||
edge.lr->computeBasis(fe.u, fe.v, spline, basis_el);
|
||||
SplineUtils::extractBasis(spline,fe.N,dNdu);
|
||||
}
|
||||
|
||||
// Assemble into matrix A
|
||||
for (size_t il=0; il<edge.MNPC[i].size(); il++) // local i-index
|
||||
if (edge.MNPC[i][il] != -1)
|
||||
{
|
||||
size_t ig = edge.MNPC[i][il]+1; // global i-index
|
||||
for (size_t jl=0; jl<edge.MNPC[i].size(); jl++) // local j-index
|
||||
if (edge.MNPC[i][jl] != -1)
|
||||
{
|
||||
size_t jg = edge.MNPC[i][jl]+1; // global j-index
|
||||
A(ig,jg) = A(ig,jg) + fe.N[il]*fe.N[jl]*fe.detJxW;
|
||||
}
|
||||
B(ig) = B(ig) + fe.N[il]*values(X)*fe.detJxW;
|
||||
} // end basis-function loop
|
||||
} // end gauss-point loop
|
||||
} // end element loop
|
||||
|
||||
#if SP_DEBUG > 2
|
||||
std::cout <<"---- Matrix A -----\n"<< A
|
||||
<<"-------------------"<< std::endl;
|
||||
std::cout <<"---- Vector B -----\n"<< B
|
||||
<<"-------------------"<< std::endl;
|
||||
|
||||
// dump mesh enumerations to file
|
||||
std::ofstream out("mesh.eps");
|
||||
lrspline->writePostscriptFunctionSpace(out);
|
||||
|
||||
std::cout <<"---- Edge-nodes (g2l-mapping) -----\n";
|
||||
int i=-1;
|
||||
for (auto d : edge.MLGN) {
|
||||
i++;
|
||||
std::cout << d.first << ": " << d.second << std::endl;
|
||||
}
|
||||
std::cout <<"-------------------"<< std::endl;
|
||||
std::cout <<"---- Element-nodes (-1 is interior element nodes) -----\n";
|
||||
for (auto d : edge.MNPC) {
|
||||
for (int c : d)
|
||||
std::cout << c << " ";
|
||||
std::cout << std::endl;
|
||||
}
|
||||
std::cout <<"-------------------"<< std::endl;
|
||||
#endif
|
||||
|
||||
// Solve the edge-global equation system
|
||||
if (!A.solve(B)) return false;
|
||||
|
||||
#if SP_DEBUG > 2
|
||||
std::cout <<"---- SOLUTION -----\n"<< B
|
||||
<<"-------------------"<< std::endl;
|
||||
#endif
|
||||
|
||||
// Store the control-point values of the projected field
|
||||
result.resize(n);
|
||||
for (size_t i = 0; i < n; i++)
|
||||
result[i] = B[i];
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool ASMu2D::edgeL2projection (const DirichletEdge& edge,
|
||||
const VecFunc& values,
|
||||
Real2DMat& result,
|
||||
double time) const
|
||||
{
|
||||
size_t n = edge.MLGN.size();
|
||||
SparseMatrix A(SparseMatrix::SUPERLU);
|
||||
std::vector<StdVector> B(nsd, StdVector(n));
|
||||
A.resize(n,n);
|
||||
|
||||
// Get Gaussian quadrature points and weights
|
||||
const double* xg = GaussQuadrature::getCoord(nGauss);
|
||||
const double* wg = GaussQuadrature::getWeight(nGauss);
|
||||
if (!xg || !wg) return false;
|
||||
|
||||
// find the normal and tangent direction for the edge
|
||||
int edgeDir, t1, t2;
|
||||
switch (edge.edg)
|
||||
{ // id normal tangent
|
||||
case LR::WEST: edgeDir = -1; t1=1; t2=2; break;
|
||||
case LR::EAST: edgeDir = +1; t1=1; t2=2; break;
|
||||
case LR::SOUTH: edgeDir = -2; t1=2; t2=1; break;
|
||||
case LR::NORTH: edgeDir = +2; t1=2; t2=1; break;
|
||||
default: return false;
|
||||
}
|
||||
|
||||
std::array<Vector,2> gpar;
|
||||
for (int d = 0; d < 2; d++)
|
||||
if (-1-d == edgeDir)
|
||||
{
|
||||
gpar[d].resize(nGauss);
|
||||
gpar[d].fill(d == 0 ? lrspline->startparam(0) : lrspline->startparam(1));
|
||||
}
|
||||
else if (1+d == edgeDir)
|
||||
{
|
||||
gpar[d].resize(nGauss);
|
||||
gpar[d].fill(d == 0 ? lrspline->endparam(0) : lrspline->endparam(1));
|
||||
}
|
||||
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
Vec3 normal;
|
||||
int iGp = 0;
|
||||
|
||||
// === Assembly loop over all elements on the patch edge =====================
|
||||
for (size_t i=0; i<edge.MLGE.size(); i++) // for all edge elements
|
||||
{
|
||||
FiniteElement fe(edge.MNPC[i].size());
|
||||
fe.iel = edge.MLGE[i];
|
||||
|
||||
// Get element edge length in the parameter space
|
||||
double dS = this->getParametricLength(fe.iel+1,t1);
|
||||
if (dS < 0.0) return false; // topology error (probably logic error)
|
||||
|
||||
// Set up control point coordinates for current element
|
||||
if (!this->getElementCoordinates(Xnod,fe.iel+1)) return false;
|
||||
|
||||
// Initialize element quantities
|
||||
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
|
||||
|
||||
// Get integration gauss points over this element
|
||||
this->getGaussPointParameters(gpar[t2-1],t2-1,nGauss,fe.iel+1,xg);
|
||||
|
||||
// --- Integration loop over all Gauss points along the edge -------------
|
||||
for (int j = 0; j < nGauss; j++)
|
||||
{
|
||||
fe.iGP = iGp++; // Global integration point counter
|
||||
|
||||
// Local element coordinates and parameter values
|
||||
// of current integration point
|
||||
fe.xi = xg[j];
|
||||
fe.eta = xg[j];
|
||||
fe.u = gpar[0][j];
|
||||
fe.v = gpar[1][j];
|
||||
|
||||
// Evaluate basis function derivatives at current integration points
|
||||
Go::BasisDerivsSf spline;
|
||||
lrspline->computeBasis(fe.u, fe.v, spline, fe.iel);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
SplineUtils::extractBasis(spline,fe.N,dNdu);
|
||||
|
||||
// Compute basis function derivatives and the edge normal
|
||||
fe.detJxW = utl::Jacobian(Jac,normal,fe.dNdX,Xnod,dNdu,t1,t2);
|
||||
if (fe.detJxW == 0.0) continue; // skip singular points
|
||||
|
||||
if (edgeDir < 0) normal *= -1.0;
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.t = time;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
fe.detJxW *= 0.5*dS*wg[j];
|
||||
|
||||
// For mixed basis, we need to compute functions separate from geometry
|
||||
if (edge.lr != lrspline.get())
|
||||
{
|
||||
// different lrspline instances enumerate elements differently
|
||||
int basis_el = edge.lr->getElementContaining(fe.u, fe.v);
|
||||
edge.lr->computeBasis(fe.u, fe.v, spline, basis_el);
|
||||
SplineUtils::extractBasis(spline,fe.N,dNdu);
|
||||
}
|
||||
|
||||
// Assemble into matrix A
|
||||
for (size_t il=0; il<edge.MNPC[i].size(); il++) // local i-index
|
||||
if (edge.MNPC[i][il] != -1)
|
||||
{
|
||||
size_t ig = edge.MNPC[i][il]+1; // global i-index
|
||||
for (size_t jl=0; jl<edge.MNPC[i].size(); jl++) // local j-index
|
||||
if (edge.MNPC[i][jl] != -1)
|
||||
{
|
||||
size_t jg = edge.MNPC[i][jl]+1; // global j-index
|
||||
A(ig,jg) = A(ig,jg) + fe.N[il]*fe.N[jl]*fe.detJxW;
|
||||
}
|
||||
Vec3 val = values(X);
|
||||
for(size_t k=0; k<nsd; k++)
|
||||
B[k](ig) = B[k](ig) + fe.N[il]*val[k]*fe.detJxW;
|
||||
} // end basis-function loop
|
||||
} // end gauss-point loop
|
||||
} // end element loop
|
||||
|
||||
#if SP_DEBUG > 2
|
||||
std::cout <<"---- Matrix A -----\n"<< A
|
||||
<<"-------------------"<< std::endl;
|
||||
for(size_t k=0; k<nsd; k++)
|
||||
{
|
||||
std::cout <<"--- Vector B(" << (k+1) << ") ----\n" << B[k]
|
||||
<<"-------------------"<< std::endl;
|
||||
}
|
||||
|
||||
// dump mesh enumerations to file
|
||||
std::ofstream out("mesh.eps");
|
||||
lrspline->writePostscriptFunctionSpace(out);
|
||||
|
||||
std::cout <<"---- Edge-nodes (g2l-mapping) -----\n";
|
||||
int i=-1;
|
||||
for (auto d : edge.MLGN) {
|
||||
i++;
|
||||
std::cout << d.first << ": " << d.second << std::endl;
|
||||
}
|
||||
std::cout <<"-------------------"<< std::endl;
|
||||
std::cout <<"---- Element-nodes (-1 is interior element nodes) -----\n";
|
||||
for (auto d : edge.MNPC) {
|
||||
for (int c : d)
|
||||
std::cout << c << " ";
|
||||
std::cout << std::endl;
|
||||
}
|
||||
std::cout <<"-------------------"<< std::endl;
|
||||
#endif
|
||||
|
||||
// Solve the edge-global equation system
|
||||
if (!A.solve(B[0], true)) return false;
|
||||
// Solve the system for the rest of the right-hand-side components (re-use LU factorization)
|
||||
for (size_t k=1; k<nsd; k++)
|
||||
if (!A.solve(B[k], false)) return false;
|
||||
|
||||
#if SP_DEBUG > 2
|
||||
for(size_t k=0; k<nsd; k++)
|
||||
{
|
||||
std::cout <<"--- Solution(" << (k+1) << ") ----\n" << B[k]
|
||||
<<"-------------------"<< std::endl;
|
||||
}
|
||||
#endif
|
||||
result.resize(nsd);
|
||||
for(size_t i=0; i<nsd; i++)
|
||||
{
|
||||
result[i].resize(n);
|
||||
std::copy(B[i].begin(), B[i].end(), result[i].begin());
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user