From 1a346dba8b058c8e64ab66fcb289c308014517aa Mon Sep 17 00:00:00 2001 From: Kjetil Andre Johannessen Date: Fri, 25 Nov 2016 14:42:42 +0100 Subject: [PATCH] Non-homogenous dirichlet BC for 2D LR-splines --- src/ASM/LR/ASMu2D.C | 200 +++++++++++++-------- src/ASM/LR/ASMu2D.h | 50 ++++-- src/ASM/LR/ASMu2Drecovery.C | 343 ++++++++++++++++++++++++++++++++++++ 3 files changed, 508 insertions(+), 85 deletions(-) diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 805c3a77..510ad9ee 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -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 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 edgeFunctions; + std::vector 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 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 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 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; ilrspline.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& func, const std::map& vfunc, double time, const std::map* g2l) { - std::map::const_iterator fit; - std::map::const_iterator vfit; + + std::map::const_iterator fit; + std::map::const_iterator vfit; + std::map::const_iterator nit; std::vector::const_iterator dit; - std::vector::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); } diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h index b996fa5a..d59c1182 100644 --- a/src/ASM/LR/ASMu2D.h +++ b/src/ASM/LR/ASMu2D.h @@ -16,6 +16,7 @@ #include "ASMunstruct.h" #include "ASM2D.h" +#include "LRSpline/LRSpline.h" #include 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 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 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 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 lrspline; //!< Pointer to the LR-spline surface object diff --git a/src/ASM/LR/ASMu2Drecovery.C b/src/ASM/LR/ASMu2Drecovery.C index 34b2f27c..8e3dda2d 100644 --- a/src/ASM/LR/ASMu2Drecovery.C +++ b/src/ASM/LR/ASMu2Drecovery.C @@ -23,8 +23,11 @@ #include "Utilities.h" #include "Profiler.h" #include "IntegrandBase.h" +#include "FiniteElement.h" #include +#include + 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 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; igetParametricLength(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 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 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 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; igetParametricLength(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 2 + std::cout <<"---- Matrix A -----\n"<< A + <<"-------------------"<< std::endl; + for(size_t k=0; kwritePostscriptFunctionSpace(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 2 + for(size_t k=0; k