From 85aaeb74ebef6e39a2483e5b68f281058d85cce3 Mon Sep 17 00:00:00 2001 From: kjetijo Date: Tue, 20 Sep 2011 09:45:59 +0000 Subject: [PATCH] added: LR splines git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1149 e10b68d5-8a6e-419e-a041-bce267b0401d --- CMakeLists.txt | 15 + cmake/Modules/FindLRSpline.cmake | 28 + src/ASM/LR/ASMu2D.C | 1365 ++++++++++++++++++++++++++++++ src/ASM/LR/ASMu2D.h | 367 ++++++++ src/ASM/LR/ASMunstruct.C | 34 + src/ASM/LR/ASMunstruct.h | 65 ++ src/SIM/SIM2D.C | 138 ++- src/SIM/SIMbase.h | 2 +- 8 files changed, 1971 insertions(+), 43 deletions(-) create mode 100644 cmake/Modules/FindLRSpline.cmake create mode 100644 src/ASM/LR/ASMu2D.C create mode 100644 src/ASM/LR/ASMu2D.h create mode 100644 src/ASM/LR/ASMunstruct.C create mode 100644 src/ASM/LR/ASMunstruct.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 4c590c5e..768827d1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,6 +51,9 @@ FIND_PACKAGE(VTFWriter) IF(NOT "${DISABLE_HDF5}") FIND_PACKAGE(HDF5) ENDIF(NOT "${DISABLE_HDF5}") +IF(NOT "${DISABLE_LRSPLINE}") + FIND_PACKAGE(LRSpline) +ENDIF(NOT "${DISABLE_LRSPLINE}") # Required libraries SET(DEPLIBS ${GoTrivariate_LIBRARIES} ${GoTools_LIBRARIES} @@ -131,6 +134,15 @@ IF(HDF5_LIBRARIES AND HDF5_INCLUDE_DIR) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHAS_HDF5=1") ENDIF(HDF5_LIBRARIES AND HDF5_INCLUDE_DIR) +IF(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS) + SET(DEPLIBS ${DEPLIBS} ${LRSpline_LIBRARIES}) + SET(INCLUDES + ${INCLUDES} + ${LRSpline_INCLUDE_DIRS} + ${PROJECT_SOURCE_DIR}/src/ASM/LR) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHAS_LRSPLINE=1") +ENDIF(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS) + INCLUDE_DIRECTORIES(${INCLUDES}) SET(EXECUTABLE_OUTPUT_PATH bin) @@ -146,6 +158,9 @@ ENDIF(NOT WIN32) # Make the IFEM library FILE(GLOB_RECURSE IFEM_SRCS ${PROJECT_SOURCE_DIR}/src/*.[Cf] ${PROJECT_SOURCE_DIR}/3rdparty/*.C) +IF(NOT(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS)) + string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/src/ASM/LR/[^;]*" "" IFEM_SRCS "${IFEM_SRCS}") +ENDIF(NOT(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS)) ADD_LIBRARY(IFEM ${IFEM_SRCS}) # Make some Apps diff --git a/cmake/Modules/FindLRSpline.cmake b/cmake/Modules/FindLRSpline.cmake new file mode 100644 index 00000000..04febd4d --- /dev/null +++ b/cmake/Modules/FindLRSpline.cmake @@ -0,0 +1,28 @@ +IF(LRSpline_INCLUDE_DIRS AND LRSpline_LIBRARIES) + SET(LRSpline_FIND_QUIETLY TRUE) +ENDIF(LRSpline_INCLUDE_DIRS AND LRSpline_LIBRARIES) + +UNSET(LRSpline_INCLUDE_DIRS CACHE) +UNSET(LRSpline_LIBRARIES CACHE) + +FIND_PATH(LRSpline_INCLUDE_DIRS + NAMES LRSpline/LRSplineSurface.h + PATHS "$ENV{HOME}/include" + "$ENV{HOME}/install/include" + /sima/libs/LRSpline/include +) + +FIND_LIBRARY(LRSpline_LIBRARIES + NAMES LRSpline + PATHS "$ENV{HOME}/lib" + "$ENV{HOME}/install/lib" + /sima/libs/LRSpline/lib + PATH_SUFFIXES LRSpline +) + +INCLUDE(FindPackageHandleStandardArgs) +IF(LRSpline_LIBRARIES) + find_package_handle_standard_args(LRSpline DEFAULT_MSG + LRSpline_INCLUDE_DIRS LRSpline_LIBRARIES) +ENDIF(LRSpline_LIBRARIES) + diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C new file mode 100644 index 00000000..ffadfe16 --- /dev/null +++ b/src/ASM/LR/ASMu2D.C @@ -0,0 +1,1365 @@ +// $Id: ASMu2D.C 1108 2011-08-23 14:51:03Z kjetijo $ +//============================================================================== +//! +//! \file ASMu2D.C +//! +//! \date September 2011 +//! +//! \author Kjetil Andre Johannessen / SINTEF +//! +//! \brief Driver for assembly of unstructured 2D spline FE models. +//! +//============================================================================== + +#include "GoTools/geometry/ObjectHeader.h" +#include "GoTools/geometry/SplineSurface.h" +#include "GoTools/geometry/SurfaceInterpolator.h" + +#include "LRSpline/LRSplineSurface.h" +#include "LRSpline/Element.h" +#include "LRSpline/Basisfunction.h" +#include + +#include "ASMu2D.h" +#include "TimeDomain.h" +#include "FiniteElement.h" +#include "GlobalIntegral.h" +#include "IntegrandBase.h" +#include "CoordinateMapping.h" +#include "GaussQuadrature.h" +#include "ElementBlock.h" +#include "Utilities.h" +#include "Profiler.h" +#include "Vec3Oper.h" +#include +#include + +// #define SP_DEBUG 5 + +ASMu2D::ASMu2D (const char* fName, unsigned char n_s, unsigned char n_f) + : ASMunstruct(2,n_s,n_f), lrspline(0) +{ + if (fName) + { + std::cout <<"\nReading patch file "<< fName << std::endl; + std::ifstream is(fName); + if (!is.good()) + std::cerr <<" *** ASMu2D: Failure opening patch file"<< std::endl; + else + this->read(is); + } +} + + +ASMu2D::ASMu2D (std::istream& is, unsigned char n_s, unsigned char n_f) + : ASMunstruct(2,n_s,n_f), lrspline(0) +{ + this->read(is); +} + + +bool ASMu2D::read (std::istream& is) +{ + if (lrspline) delete lrspline; + + Go::ObjectHeader head; + Go::SplineSurface *surf = new Go::SplineSurface(); + is >> head >> *surf; + + lrspline = new LR::LRSplineSurface(surf); + + // Eat white-space characters to see if there is more data to read + char c; + while (is.get(c)) + if (!isspace(c)) + { + is.putback(c); + break; + } + + if (!is.good() && !is.eof()) + { + std::cerr <<" *** ASMu2D::read: Failure reading spline data"<< std::endl; + delete lrspline; + lrspline = 0; + return false; + } + else if (lrspline->dimension() < 2) + { + std::cerr <<" *** ASMu2D::read: Invalid spline lrsplineace patch, dim=" + << lrspline->dimension() << std::endl; + delete lrspline; + lrspline = 0; + return false; + } + else if (lrspline->dimension() < nsd) + { + std::cout <<" ** ASMu2D::read: The dimension of this lrsplineace patch " + << lrspline->dimension() <<" is less than nsd="<< nsd + <<".\n Resetting nsd to "<< lrspline->dimension() + <<" for this patch."<< std::endl; + nsd = lrspline->dimension(); + } + + geo = lrspline; + return true; +} + + +bool ASMu2D::write (std::ostream& os, int) const +{ + if (!lrspline) return false; + + os <<"200 1 0 0\n"; + os << *lrspline; + + return os.good(); +} + + +void ASMu2D::clear () +{ + // Erase spline data + if (lrspline) delete lrspline; + lrspline = 0; + geo = 0; + + // Erase the FE data + ASMbase::clear(); +} + +bool ASMu2D::cornerRefine (int minBasisfunctions) +{ + if (!lrspline ) return false; + + double h = 1.0; + int nBasis = lrspline->nBasisFunctions(); + double unif_step_h = 1.0 / ((minBasisfunctions - nBasis) / 3.0 + 1.0); + while(lrspline->nBasisFunctions() < minBasisfunctions) { + lrspline->insert_const_u_edge(h-unif_step_h, 0, h); + lrspline->insert_const_v_edge(h-unif_step_h, 0, h); + h -= unif_step_h; + } + std::ofstream meshFile; + meshFile.open("mesh.eps"); + lrspline->writePostscriptMesh(meshFile); + meshFile.close(); + return true; +} + +bool ASMu2D::diagonalRefine (int minBasisfunctions) +{ + if (!lrspline ) return false; + + double end1 = lrspline->endparam_u(); + double end2 = lrspline->endparam_v(); + double h = 1.0; + int iter = 0; + double u = h/2.0; + double v = h/2.0; + while(lrspline->nBasisFunctions() < minBasisfunctions) { + lrspline->insert_const_u_edge(u, (iter-1<0) ? 0 : (iter-1)*h, ((iter+2)*h>end2) ? end2 : (iter+2)*h); + lrspline->insert_const_v_edge(v, (iter-1<0) ? 0 : (iter-1)*h, ((iter+2)*h>end1) ? end1 : (iter+2)*h); + u += h; + v += h; + iter++; + if( u>end1 ) { + h /= 2.0; + iter = 0; + u = h/2.0; + v = h/2.0; + } + } + std::ofstream meshFile; + meshFile.open("mesh.eps"); + lrspline->writePostscriptMesh(meshFile); + meshFile.close(); + return true; +} + +bool ASMu2D::uniformRefine (int minBasisfunctions) +{ + if (!lrspline ) return false; + + double end1 = lrspline->endparam_u(); + double end2 = lrspline->endparam_v(); + double h = 1.0; + bool step_u = true; + double u = h/2.0; + double v = h/2.0; + while(lrspline->nBasisFunctions() < minBasisfunctions) { + if(step_u) { + lrspline->insert_const_u_edge(u, 0, end2); + u += h; + if(u > end1) { + step_u = !step_u; + u = h/4.0; + } + } else { + lrspline->insert_const_v_edge(v, 0, end1); + v += h; + if(v > end2) { + step_u = !step_u; + v = h/4.0; + h /= 2.0; + } + } + } + std::ofstream meshFile; + meshFile.open("mesh.eps"); + lrspline->writePostscriptMesh(meshFile); + meshFile.close(); + return true; +} + + +bool ASMu2D::generateFEMTopology () +{ + if (!lrspline) return false; + + const int nBasis = lrspline->nBasisFunctions(); + const int nElements = lrspline->nElements(); + + const int p1 = lrspline->order_u(); + const int p2 = lrspline->order_v(); +#ifdef SP_DEBUG + lrspline->write(std::cout); +#endif + // Consistency checks, just to be fool-proof + if (nBasis < 4) return false; + if (p1 < 1 || p1 < 1) return false; + if (p1*p2 > nBasis) return false; + + MLGE.resize(nElements,0); + MLGN.resize(nBasis); + MNPC.resize(MLGE.size()); + lrspline->generateIDs(); + + int iel = 0; + std::vector::iterator el_it = lrspline->elementBegin(); + std::vector::iterator b_it; + for(iel=0; ielnBasisFunctions(); + MLGE[iel] = ++gEl; // global element number over all patches + MNPC[iel].resize(nSupportFunctions, 0); + + b_it = (*el_it)->supportBegin(); + for(int lnod=0; lnodgetId(); + } + + for(size_t i=0; i= nnodI-1) + { + indxI = 1; + iinod += inc[1] - inc[0]*(nnodI-2); + } + + return ret; +} + + +bool ASMu2D::assignNodeNumbers () +{ + for(size_t i=0; iconnectBasis(edge,neighbor,nedge,revers); +} + + +bool ASMu2D::connectBasis (int edge, ASMu2D& neighbor, int nedge, bool revers, + int basis, int slave, int master) +{ + // Set up the slave node numbers for this surface patch + + int n1, n2; + if (!this->getSize(n1,n2,basis)) return false; + int node = slave+1, i1 = 0; + + switch (edge) + { + case 2: // Positive I-direction + node += n1-1; + case 1: // Negative I-direction + i1 = n1; + n1 = n2; + break; + + case 4: // Positive J-direction + node += n1*(n2-1); + case 3: // Negative J-direction + i1 = 1; + break; + + default: + std::cerr <<" *** ASMu2D::connectPatch: Invalid slave edge " + << edge << std::endl; + return false; + } + + int i; + IntVec slaveNodes(n1,0); + for (i = 0; i < n1; i++, node += i1) + slaveNodes[i] = node; + + // Set up the master node numbers for the neighboring surface patch + + if (!neighbor.getSize(n1,n2,basis)) return false; + node = master+1; i1 = 0; + + switch (nedge) + { + case 2: // Positive I-direction + node += n1-1; + case 1: // Negative I-direction + i1 = n1; + n1 = n2; + break; + + case 4: // Positive J-direction + node += n1*(n2-1); + case 3: // Negative J-direction + i1 = 1; + break; + + default: + std::cerr <<" *** ASMu2D::connectPatch: Invalid master edge " + << nedge << std::endl; + return false; + } + + if (n1 != (int)slaveNodes.size()) + { + std::cerr <<" *** ASMu2D::connectPatch: Non-matching edges, sizes " + << n1 <<" and "<< slaveNodes.size() << std::endl; + return false; + } + + const double xtol = 1.0e-4; + for (i = 0; i < n1; i++, node += i1) + { + int k = revers ? n1-i-1 : i; + if (!neighbor.getCoord(node).equal(this->getCoord(slaveNodes[k]),xtol)) + { + std::cerr <<" *** ASMu2D::connectPatch: Non-matching nodes " + << node <<": "<< neighbor.getCoord(node) + <<"\n and " + << slaveNodes[k] <<": "<< this->getCoord(slaveNodes[k]) + << std::endl; + return false; + } + else + ASMbase::collapseNodes(neighbor.MLGN[node-1],MLGN[slaveNodes[k]-1]); + } + + return true; +} + + +void ASMu2D::closeEdges (int dir, int basis, int master) +{ + int n1, n2; + if (basis < 1) basis = 1; + if (!this->getSize(n1,n2,basis)) return; + + switch (dir) + { + case 1: // Edges are closed in I-direction + for (int i2 = 1; i2 <= n2; i2++, master += n1) + this->makePeriodic(master,master+n1-1); + break; + + case 2: // Edges are closed in J-direction + for (int i1 = 1; i1 <= n1; i1++, master++) + this->makePeriodic(master,master+n1*(n2-1)); + break; + } +} + +#endif + + +void ASMu2D::constrainEdge (int dir, int dof, int code) +{ + std::vector edgeFunctions; + switch (dir) + { + case 1: // Right edge (positive I-direction) + lrspline->getEdgeFunctions(edgeFunctions, LR::EAST); + break; + case -1: // Left edge (negative I-direction) + lrspline->getEdgeFunctions(edgeFunctions, LR::WEST); + break; + case 2: // Back edge (positive J-direction) + lrspline->getEdgeFunctions(edgeFunctions, LR::NORTH); + break; + case -2: // Front edge (negative J-direction) + lrspline->getEdgeFunctions(edgeFunctions, LR::SOUTH); + break; + } + for(size_t i=0; iprescribe(edgeFunctions[i]->getId()+1,dof,code); +} + + +void ASMu2D::constrainCorner (int I, int J, int dof, int code) +{ + std::vector edgeFunctions; + if (I == 0 && J == 0) + lrspline->getEdgeFunctions(edgeFunctions, LR::SOUTH_WEST); + else if(I > 0 && J == 0) + lrspline->getEdgeFunctions(edgeFunctions, LR::SOUTH_EAST); + else if(I == 0 && J > 0) + lrspline->getEdgeFunctions(edgeFunctions, LR::NORTH_WEST); + else + lrspline->getEdgeFunctions(edgeFunctions, LR::NORTH_EAST); + +#ifdef SP_DEBUG + if(edgeFunctions.size() != 1) { + std::cerr <<" *** ASMu2D::constrainCorner: more than one corner" + <<" returned from lrspline->getEdgeFunctions()" << std::endl; + return; + } +#endif + + this->prescribe(edgeFunctions[0]->getId()+1,dof,code); +} + + +// Hopefully we don't have to constrain non-corner singlenodes inside patches +#if 0 +void ASMu2D::constrainNode (double xi, double eta, int dof, int code) +{ + if (xi < 0.0 || xi > 1.0) return; + if (eta < 0.0 || eta > 1.0) return; + + int n1, n2; + if (!this->getSize(n1,n2,1)) return; + + int node = 1; + if (xi > 0.0) node += int(0.5+(n1-1)*xi); + if (eta > 0.0) node += n1*int(0.5+(n2-1)*eta); + + this->prescribe(node,dof,code); +} +#endif + + +#define DERR -999.99 + +double ASMu2D::getParametricArea (int iel) const +{ +#ifdef INDEX_CHECK + if (iel < 1 || (size_t)iel > MNPC.size()) + { + std::cerr <<" *** ASMu2D::getParametricArea: Element index "<< iel + <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; + return DERR; + } +#endif + if (MNPC[iel-1].empty()) + return 0.0; + + return lrspline->getElement(iel-1)->area(); +} + + +double ASMu2D::getParametricLength (int iel, int dir) const +{ +#ifdef INDEX_CHECK + if (iel < 1 || (size_t)iel > MNPC.size()) + { + std::cerr <<" *** ASMu2D::getParametricLength: Element index "<< iel + <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; + return DERR; + } +#endif + if (MNPC[iel-1].empty()) + return 0.0; + + LR::Element *el = lrspline->getElement(iel-1); + switch (dir) + { + case 1: return el->vmax() - el->vmin(); + case 2: return el->umax() - el->umin(); + } + + std::cerr <<" *** ASMu2D::getParametricLength: Invalid edge direction " + << dir << std::endl; + return DERR; +} + + +bool ASMu2D::getElementCoordinates (Matrix& X, int iel) const +{ +#ifdef INDEX_CHECK + if (iel < 1 || (size_t)iel > MNPC.size()) + { + std::cerr <<" *** ASMu2D::getElementCoordinates: Element index "<< iel + <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; + return false; + } +#endif + + const IntVec& mnpc = MNPC[iel-1]; + X.resize(nsd,mnpc.size()); + + std::vector::iterator bit = lrspline->getElement(iel-1)->supportBegin(); + for (size_t n = 0; n < mnpc.size(); n++, bit++) + for (size_t i = 0; i < nsd; i++) + X(i+1,n+1) = (*bit)->controlpoint_[i]; + +#if SP_DEBUG > 2 + std::cout <<"\nCoordinates for element "<< iel << X << std::endl; +#endif + return true; +} + + +void ASMu2D::getNodalCoordinates (Matrix& X) const +{ + const int nBasis = lrspline->nBasisFunctions(); + X.resize(nsd,nBasis); + + std::vector::iterator bit = lrspline->basisBegin(); + + for(int inod=0; inodcontrolpoint_[i]; +} + + +Vec3 ASMu2D::getCoord (size_t inod) const +{ + Vec3 X; + LR::Basisfunction* basis = lrspline->getBasisfunction(inod-1); + for (unsigned char i = 0; i < nsd; i++) + X[i] = basis->controlpoint_[i]; + + return X; +} + + +void ASMu2D::getGaussPointParameters (Vector& uGP, int dir, int nGauss, + int iel, const double* xi) const +{ +#ifdef INDEX_CHECK + if (iel < 1 || (size_t)iel > MNPC.size()) + { + std::cerr <<" *** ASMu2D::getGaussPointParameters: Element index "<< iel + <<" out of range [1,"<< MNPC.size() <<"]."<< std::endl; + return; + } +#endif + LR::Element* el = lrspline->getElement(iel-1); + + uGP.resize(nGauss); + double ustart = (dir==0) ? el->umin() : el->vmin(); + double ustop = (dir==0) ? el->umax() : el->vmax(); + + for (int i = 1; i <= nGauss; i++) + uGP(i) = 0.5*((ustop-ustart)*xi[i-1] + ustop+ustart); +} + + +/*! + \brief Computes the characteristic element length from nodal coordinates. +*/ + +#if 0 +static double getElmSize (int p1, int p2, const Matrix& X) +{ + int n = X.rows(); + int i, j, id1, id2; + double value, v1, h = 1.0e12; + + // Y-direction + for (i = 1; i <= p1; i++) + { + id1 = i; + id2 = id1 + (p2-1)*p1; + value = 0.0; + for (j = 1; j <= n; j++) + { + v1 = X(j,id2) - X(j,id1); + value += v1*v1; + } + if (value < h) h = value; + } + + // X-direction + for (j = 0; j < p2; j++) + { + id1 = j*p1 + 1; + id2 = id1 + p1 - 1; + value = 0.0; + for (i = 1; i <= n; i++) + { + v1 = X(i,id2) - X(i,id1); + value += v1*v1; + } + if (value < h) h = value; + } + + return sqrt(h); +} +#endif + + +void ASMu2D::extractBasis (const Go::BasisDerivsSf& spline, + Vector& N, Matrix& dNdu) +{ + dNdu.resize(N.size(),2); + + size_t jp, n = 1; + for (jp = 0; jp < N.size(); jp++, n++) + { + N (n) = spline.basisValues[jp]; + dNdu(n,1) = spline.basisDerivs_u[jp]; + dNdu(n,2) = spline.basisDerivs_v[jp]; + } +} + + +void ASMu2D::extractBasis (const Go::BasisDerivsSf2& spline, + Vector& N, Matrix& dNdu, Matrix3D& d2Ndu2) +{ + dNdu .resize(N.size(),2); + d2Ndu2.resize(N.size(),2,2); + + size_t jp, n = 1; + for (jp = 0; jp < N.size(); jp++, n++) + { + N (n) = spline.basisValues[jp]; + dNdu (n,1) = spline.basisDerivs_u[jp]; + dNdu (n,2) = spline.basisDerivs_v[jp]; + d2Ndu2(n,1,1) = spline.basisDerivs_uu[jp]; + d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline.basisDerivs_uv[jp]; + d2Ndu2(n,2,2) = spline.basisDerivs_vv[jp]; + } +} + + +bool ASMu2D::integrate (Integrand& integrand, + GlobalIntegral& glInt, + const TimeDomain& time, + const LintegralVec& locInt) +{ + if (!lrspline) return true; // silently ignore empty patches + + PROFILE2("ASMu2D::integrate(I)"); + + // Get Gaussian quadrature points and weights + const double* xg = GaussQuadrature::getCoord(nGauss); + const double* wg = GaussQuadrature::getWeight(nGauss); + if (!xg || !wg) return false; + + // Get the reduced integration quadrature points, if needed + const double* xr = 0; + int nRed = integrand.getIntegrandType() - 10; + if (nRed < 1) + nRed = nRed < 0 ? nGauss : 0; + else if (!(xr = GaussQuadrature::getCoord(nRed))) + return false; + + Matrix dNdu, Xnod, Jac; + Matrix3D d2Ndu2, Hess; + Vec4 X; + + + // === Assembly loop over all elements in the patch ========================== + + std::vector::iterator el; + int iel = 1; + for(el = lrspline->elementBegin(); el!=lrspline->elementEnd(); el++, iel++) + { + FiniteElement fe(MNPC[iel-1].size()); + fe.iel = MLGE[iel-1]; + if (fe.iel < 1) continue; // zero-area element + + // Get element area in the parameter space + double dA = this->getParametricArea(iel); + if (dA < 0.0) return false; // topology error (probably logic error) + + // Set up control point (nodal) coordinates for current element + if (!this->getElementCoordinates(Xnod,iel)) return false; + + // Compute parameter values of the Gauss points over this element + Vector gpar[2], redpar[2]; + for (int d = 0; d < 2; d++) + { + this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg); + if (integrand.getIntegrandType() > 10) + this->getGaussPointParameters(redpar[d],d,nRed,iel,xr); + } + + +#if 0 + // Compute characteristic element length, if needed + if (integrand.getIntegrandType() == 2) + fe.h = getElmSize(p1,p2,Xnod); + + else if (integrand.getIntegrandType() == 3) + { + // --- Compute average value of basis functions over the element ------- + + fe.Navg.resize(p1*p2,true); + double area = 0.0; + int ip = ((i2-p2)*nGauss*nel1 + i1-p1)*nGauss; + for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1)) + for (int i = 0; i < nGauss; i++, ip++) + { + // Fetch basis function derivatives at current integration point + extractBasis(spline[ip],fe.N,dNdu); + + // Compute Jacobian determinant of coordinate mapping + // and multiply by weight of current integration point + double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false); + double weight = 0.25*dA*wg[i]*wg[j]; + + // Numerical quadrature + fe.Navg.add(fe.N,detJac*weight); + area += detJac*weight; + } + + // Divide by element area + fe.Navg /= area; + } + + else if (integrand.getIntegrandType() == 4) + { + // Compute the element center + Go::Point X0; + double u0 = 0.5*(gpar[0](1,i1-p1+1) + gpar[0](nGauss,i1-p1+1)); + double v0 = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](nGauss,i2-p2+1)); + lrspline->point(X0,u0,v0); + for (unsigned char i = 0; i < nsd; i++) + X[i] = X0[i]; + } +#endif + + // Initialize element quantities + if (!integrand.initElement(MNPC[iel-1],X,nRed*nRed)) return false; + + // Caution: Unless locInt is empty, we assume it points to an array of + // LocalIntegral pointers, of length at least the number of elements in + // the model (as defined by the highest number in the MLGE array). + // If the array is shorter than this, expect a segmentation fault. + LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1]; + + +#if 0 + if (integrand.getIntegrandType() > 10) + { + // --- Selective reduced integration loop ------------------------------ + + int ip = ((i2-p2)*nRed*nel1 + i1-p1)*nRed; + for (int j = 0; j < nRed; j++, ip += nRed*(nel1-1)) + for (int i = 0; i < nRed; i++, ip++) + { + // Local element coordinates of current integration point + fe.xi = xr[i]; + fe.eta = xr[j]; + + // Parameter values of current integration point + fe.u = redpar[0](i+1,i1-p1+1); + fe.v = redpar[1](j+1,i2-p2+1); + + // Fetch basis function derivatives at current point + extractBasis(splineRed[ip],fe.N,dNdu); + + // Compute Jacobian inverse and derivatives + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + + // Compute the reduced integration terms of the integrand + if (!integrand.reducedInt(fe)) + return false; + } + } +#endif + + + // --- Integration loop over all Gauss points in each direction ---------- + + for (int j = 0; j < nGauss; j++) + for (int i = 0; i < nGauss; i++) + { + // Local element coordinates of current integration point + fe.xi = xg[i]; + fe.eta = xg[j]; + + // Parameter values of current integration point + fe.u = gpar[0](i+1); + fe.v = gpar[1](j+1); + + // compute basis function values at integration point + Go::BasisDerivsSf spline; + Go::BasisDerivsSf2 spline2; + Go::BasisDerivsSf splineRed; +// I'll impement second order differentiation in LRSplines soon - promise +#if 0 + if (integrand.getIntegrandType() == 2) + lrspline->computeBasis(fe.u,fe.v,spline2, iel); + else +#endif + lrspline->computeBasis(fe.u,fe.v,spline, iel-1); + if (integrand.getIntegrandType() > 10) + lrspline->computeBasis(fe.u,fe.v,splineRed, iel-1); + + // Fetch basis function derivatives at current integration point + if (integrand.getIntegrandType() == 2) + extractBasis(spline2,fe.N,dNdu,d2Ndu2); + else + extractBasis(spline,fe.N,dNdu); + +#if SP_DEBUG > 4 + std::cout <<"\nBasis functions at a integration point "; + std::cout <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1] + <<" left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1] << std::endl; + for (size_t ii = 0; ii < spline.basisValues.size(); ii++) + std::cout << 1+ii <<'\t'<< spline.basisValues[ii] <<'\t' + << spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii] << std::endl; +#endif + + + // Compute Jacobian inverse of coordinate mapping and derivatives + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + if (fe.detJxW == 0.0) continue; // skip singular points + + // Compute Hessian of coordinate mapping and 2nd order derivatives + if (integrand.getIntegrandType() == 2) + if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) + return false; + +#if SP_DEBUG > 4 + std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl; +#endif + + // Cartesian coordinates of current integration point + X = Xnod * fe.N; + X.t = time.t; + + // Evaluate the integrand and accumulate element contributions + fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + if (!integrand.evalInt(elmInt,fe,time,X)) + return false; + } + + // Finalize the element quantities + if (!integrand.finalizeElement(elmInt,time)) + return false; + + // Assembly of global system integral + if (!glInt.assemble(elmInt,fe.iel)) + return false; + } + + return true; +} + + +bool ASMu2D::integrate (Integrand& integrand, int lIndex, + GlobalIntegral& glInt, + const TimeDomain& time, + const LintegralVec& locInt) +{ + if (!lrspline) return true; // silently ignore empty patches + +#if 0 + PROFILE2("ASMu2D::integrate(B)"); + + // 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 parametric direction of the edge normal {-2,-1, 1, 2} + const int edgeDir = (lIndex+1)/(lIndex%2 ? -2 : 2); + + const int t1 = abs(edgeDir); // Tangent direction normal to the patch edge + const int t2 = 3-abs(edgeDir); // Tangent direction along the patch edge + + // Compute parameter values of the Gauss points along the whole patch edge + Matrix gpar[2]; + for (int d = 0; d < 2; d++) + if (-1-d == edgeDir) + { + gpar[d].resize(1,1); + gpar[d].fill(d == 0 ? lrspline->startparam_u() : lrspline->startparam_v()); + } + else if (1+d == edgeDir) + { + gpar[d].resize(1,1); + gpar[d].fill(d == 0 ? lrspline->endparam_u() : lrspline->endparam_v()); + } + else + this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg); + + // Evaluate basis function derivatives at all integration points + std::vector spline; + lrspline->computeBasisGrid(gpar[0],gpar[1],spline); + + const int p1 = lrspline->order_u(); + const int p2 = lrspline->order_v(); + const int n1 = lrspline->numCoefs_u(); + const int n2 = lrspline->numCoefs_v(); + + FiniteElement fe(p1*p2); + fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; + fe.u = gpar[0](1,1); + fe.v = gpar[1](1,1); + + Matrix dNdu, Xnod, Jac; + Vec4 X; + Vec3 normal; + + + // === Assembly loop over all elements on the patch edge ===================== + + std::vector::iterator el + int iel = 1; + for(el = lrspline->elementBegin(); el!=lrspline->elementEnd(); el++, iel++) + { + fe.iel = MLGE[iel-1]; + if (fe.iel < 1) continue; // zero-area element + + // Skip elements that are not on current boundary edge + bool skipMe = false; + switch (edgeDir) + { + case -1: if (i1 > p1) skipMe = true; break; + case 1: if (i1 < n1) skipMe = true; break; + case -2: if (i2 > p2) skipMe = true; break; + case 2: if (i2 < n2) skipMe = true; break; + } + if (skipMe) continue; + + // Get element edge length in the parameter space + double dS = this->getParametricLength(iel,t1); + if (dS < 0.0) return false; // topology error (probably logic error) + + // Set up control point coordinates for current element + if (!this->getElementCoordinates(Xnod,iel)) return false; + + // Initialize element quantities + if (!integrand.initElementBou(MNPC[iel-1])) return false; + + // Caution: Unless locInt is empty, we assume it points to an array of + // LocalIntegral pointers, of length at least the number of elements in + // the model (as defined by the highest number in the MLGE array). + // If the array is shorter than this, expect a segmentation fault. + LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1]; + + + // --- Integration loop over all Gauss points along the edge ------------- + + int ip = (t1 == 1 ? i2-p2 : i1-p1)*nGauss; + for (int i = 0; i < nGauss; i++, ip++) + { + // Local element coordinates and parameter values + // of current integration point + if (gpar[0].size() > 1) + { + fe.xi = xg[i]; + fe.u = gpar[0](i+1,i1-p1+1); + } + if (gpar[1].size() > 1) + { + fe.eta = xg[i]; + fe.v = gpar[1](i+1,i2-p2+1); + } + + // Fetch basis function derivatives at current integration point + extractBasis(spline[ip],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.t; + + // Evaluate the integrand and accumulate element contributions + fe.detJxW *= 0.5*dS*wg[i]; + if (!integrand.evalBou(elmInt,fe,time,X,normal)) + return false; + } + + // Assembly of global system integral + if (!glInt.assemble(elmInt,fe.iel)) + return false; + } + +#endif + return true; +} + + +int ASMu2D::evalPoint (const double* xi, double* param, Vec3& X) const +{ + if (!lrspline) return -2; + + param[0] = (1.0-xi[0])*lrspline->startparam_u() + xi[0]*lrspline->endparam_u(); + param[1] = (1.0-xi[1])*lrspline->startparam_v() + xi[1]*lrspline->endparam_v(); + + Go::Point X0; + lrspline->point(X0,param[0],param[1]); + for (unsigned char d = 0; d < nsd; d++) + X[d] = X0[d]; + + return 0; +} + + +#if 0 +bool ASMu2D::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) const +{ + if (!lrspline) return false; + + if (nSegPerSpan < 1) + { + std::cerr <<" *** ASMu2D::getGridParameters: Too few knot-span points " + << nSegPerSpan+1 <<" in direction "<< dir << std::endl; + return false; + } + + RealArray::const_iterator uit = lrspline->basis(dir).begin(); + double ucurr = 0.0, uprev = *(uit++); + while (uit != lrspline->basis(dir).end()) + { + ucurr = *(uit++); + if (ucurr > uprev) + if (nSegPerSpan == 1) + prm.push_back(uprev); + else for (int i = 0; i < nSegPerSpan; i++) + { + double xg = (double)(2*i-nSegPerSpan)/(double)nSegPerSpan; + prm.push_back(0.5*(ucurr*(1.0+xg) + uprev*(1.0-xg))); + } + uprev = ucurr; + } + + if (ucurr > prm.back()) + prm.push_back(ucurr); + return true; +} + + +bool ASMu2D::getGrevilleParameters (RealArray& prm, int dir) const +{ + if (!lrspline) return false; + + const Go::BsplineBasis& basis = lrspline->basis(dir); + + prm.resize(basis.numCoefs()); + for (size_t i = 0; i < prm.size(); i++) + prm[i] = basis.grevilleParameter(i); + + return true; +} + + +bool ASMu2D::tesselate (ElementBlock& grid, const int* npe) const +{ + // Compute parameter values of the nodal points + RealArray gpar[2]; + for (int dir = 0; dir < 2; dir++) + if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1)) + return false; + + // Evaluate the spline lrsplineace at all points + size_t nx = gpar[0].size(); + size_t ny = gpar[1].size(); + RealArray XYZ(lrspline->dimension()*nx*ny); + lrspline->gridEvaluator(XYZ,gpar[0],gpar[1]); + + // Establish the block grid coordinates + size_t i, j, l; + grid.resize(nx,ny); + for (i = l = 0; i < grid.getNoNodes(); i++, l += lrspline->dimension()) + for (j = 0; j < nsd; j++) + grid.setCoor(i,j,XYZ[l+j]); + + // Establish the block grid topology + int n[4], ip = 0; + for (j = 1, n[1] = 0; j < ny; j++) + { + n[0] = n[1]; + n[1] = n[0] + 1; + n[2] = n[1] + nx; + n[3] = n[1] + nx-1; + for (i = 1; i < nx; i++) + for (l = 0; l < 4; l++) + grid.setNode(ip++,n[l]++); + } + + return true; +} + + +void ASMu2D::scatterInd (int n1, int n2, int p1, int p2, + const int* start, IntVec& index) +{ + index.reserve(p1*p2); + int ip = ((start[1]-p2+1))*n1 + (start[0]-p1+1); + for (int i2 = 0; i2 < p2; i2++, ip += n1-p1) + for (int i1 = 0; i1 < p1; i1++, ip++) + index.push_back(ip); +} + + +bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol, + const int* npe) const +{ + // Compute parameter values of the result sampling points + RealArray gpar[2]; + for (int dir = 0; dir < 2; dir++) + if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1)) + return false; + + // Evaluate the primary solution at all sampling points + return this->evalSolution(sField,locSol,gpar); +} + + +bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol, + const RealArray* gpar, bool regular) const +{ + // Evaluate the basis functions at all points + std::vector spline; + if (regular) + lrspline->computeBasisGrid(gpar[0],gpar[1],spline); + else if (gpar[0].size() == gpar[1].size()) + { + spline.resize(gpar[0].size()); + for (size_t i = 0; i < spline.size(); i++) + lrspline->computeBasis(gpar[0][i],gpar[1][i],spline[i]); + } + else + return false; + + const int p1 = lrspline->order_u(); + const int p2 = lrspline->order_v(); + const int n1 = lrspline->numCoefs_u(); + const int n2 = lrspline->numCoefs_v(); + size_t nComp = locSol.size() / this->getNoNodes(); + if (nComp*this->getNoNodes() != locSol.size()) + return false; + + Matrix Xtmp; + Vector Ytmp; + + // Evaluate the primary solution field at each point + size_t nPoints = spline.size(); + sField.resize(nComp,nPoints); + for (size_t i = 0; i < nPoints; i++) + { + IntVec ip; + scatterInd(n1,n2,p1,p2,spline[i].left_idx,ip); + + utl::gather(ip,nComp,locSol,Xtmp); + Xtmp.multiply(spline[i].basisValues,Ytmp); + sField.fillColumn(1+i,Ytmp); + } + + return true; +} + + +bool ASMu2D::evalSolution (Matrix& sField, const Integrand& integrand, + const int* npe, bool project) const +{ + if (npe) + { + // Compute parameter values of the result sampling points + RealArray gpar[2]; + if (this->getGridParameters(gpar[0],0,npe[0]-1) && + this->getGridParameters(gpar[1],1,npe[1]-1)) + if (project) + { + // Project the secondary solution onto the spline basis + Go::SplineSurface* s = this->projectSolution(integrand); + if (s) + { + // Evaluate the projected field at the result sampling points + const Vector& svec = sField; // using utl::matrix cast operator + sField.resize(s->dimension(),gpar[0].size()*gpar[1].size()); + s->gridEvaluator(const_cast(svec),gpar[0],gpar[1]); + delete s; + return true; + } + } + else + // Evaluate the secondary solution directly at all sampling points + return this->evalSolution(sField,integrand,gpar); + } + else + { + // Project the secondary solution onto the spline basis + Go::SplineSurface* s = this->projectSolution(integrand); + if (s) + { + // Extract control point values from the spline object + sField.resize(s->dimension(),s->numCoefs_u()*s->numCoefs_v()); + sField.fill(&(*s->coefs_begin())); + delete s; + return true; + } + } + + std::cerr <<" *** ASMu2D::evalSolution: Failure!"<< std::endl; + return false; +} + + +Go::GeomObject* ASMu2D::evalSolution (const Integrand& integrand) const +{ + return this->projectSolution(integrand); +} + + +Go::SplineSurface* ASMu2D::projectSolution (const Integrand& integrand) const +{ + // Compute parameter values of the result sampling points (Greville points) + RealArray gpar[2]; + for (int dir = 0; dir < 2; dir++) + if (!this->getGrevilleParameters(gpar[dir],dir)) + return 0; + + // Evaluate the secondary solution at all sampling points + Matrix sValues; + if (!this->evalSolution(sValues,integrand,gpar)) + return 0; + + // Project the results onto the spline basis to find control point + // values based on the result values evaluated at the Greville points. + // Note that we here implicitly assume the the number of Greville points + // equals the number of control points such that we don't have to resize + // the result array. Think that is always the case, but beware if trying + // other projection schemes later. + + RealArray weights; + if (lrspline->rational()) + lrspline->getWeights(weights); + + const Vector& vec = sValues; + return Go::SurfaceInterpolator::regularInterpolation(lrspline->basis(0), + lrspline->basis(1), + gpar[0], gpar[1], + const_cast(vec), + sValues.rows(), + lrspline->rational(), + weights); +} + + +bool ASMu2D::evalSolution (Matrix& sField, const Integrand& integrand, + const RealArray* gpar, bool regular) const +{ + sField.resize(0,0); + + // Evaluate the basis functions and their derivatives at all points + std::vector spline; + if (regular) + lrspline->computeBasisGrid(gpar[0],gpar[1],spline); + else if (gpar[0].size() == gpar[1].size()) + { + spline.resize(gpar[0].size()); + std::vector tmpSpline(1); + for (size_t i = 0; i < spline.size(); i++) + { + lrspline->computeBasisGrid(RealArray(1,gpar[0][i]), + RealArray(1,gpar[1][i]), + tmpSpline); + spline[i] = tmpSpline.front(); + } + // TODO: Request a GoTools method replacing the above: + // void SplineSurface::computeBasisGrid(double param_u, double param_v, + // BasisDerivsSf& result) const + /* + for (size_t i = 0; i < spline.size(); i++) + lrspline->computeBasis(gpar[0][i],gpar[1][i],spline[i]); + */ + } + else + return false; + + const int p1 = lrspline->order_u(); + const int p2 = lrspline->order_v(); + const int n1 = lrspline->numCoefs_u(); + const int n2 = lrspline->numCoefs_v(); + + // Fetch nodal (control point) coordinates + Matrix Xnod, Xtmp; + this->getNodalCoordinates(Xnod); + + Vector N(p1*p2), solPt; + Matrix dNdu, dNdX, Jac; + + // Evaluate the secondary solution field at each point + size_t nPoints = spline.size(); + for (size_t i = 0; i < nPoints; i++) + { + // Fetch indices of the non-zero basis functions at this point + IntVec ip; + scatterInd(n1,n2,p1,p2,spline[i].left_idx,ip); + + // Fetch associated control point coordinates + utl::gather(ip,nsd,Xnod,Xtmp); + + // Fetch basis function derivatives at current integration point + extractBasis(spline[i],N,dNdu); + + // Compute the Jacobian inverse + if (utl::Jacobian(Jac,dNdX,Xtmp,dNdu) == 0.0) // Jac = (Xtmp * dNdu)^-1 + continue; // skip singular points + + // Now evaluate the solution field + if (!integrand.evalSol(solPt,N,dNdX,Xtmp*N,ip)) + return false; + else if (sField.empty()) + sField.resize(solPt.size(),nPoints,true); + + sField.fillColumn(1+i,solPt); + } + + return true; +} + +#endif + diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h new file mode 100644 index 00000000..91dcee21 --- /dev/null +++ b/src/ASM/LR/ASMu2D.h @@ -0,0 +1,367 @@ +// $Id: ASMu2D.h 1108 2011-08-23 14:51:03Z kjetijo $ +//============================================================================== +//! +//! \file ASMu2D.h +//! +//! \date September 2011 +//! +//! \author Kjetil Andre Johannessen / SINTEF +//! +//! \brief Driver for assembly of unstructured 2D spline FE models. +//! +//============================================================================== + +#ifndef _ASM_U2D_H +#define _ASM_U2D_H + +#include "ASMunstruct.h" + +namespace Go { + class SplineSurface; + class BasisDerivsSf; + class BasisDerivsSf2; +} +namespace LR { + class LRSplineSurface; +} + + +/*! + \brief Driver for assembly of structured 2D spline FE models. + \details This class contains methods common for structured 2D spline patches. +*/ + +class ASMu2D : public ASMunstruct +{ + //! \brief Struct for nodal point data. + struct IJ + { + int I; //!< Index in first parameter direction + int J; //!< Index in second parameter direction + }; + +// topology-stuff like connecting multiple patches. Will be introduced later + +#if 0 + + //! \brief Struct for edge node definitions. + struct Edge + { + int icnod; //!< Global node number of first interior point along the edge + int incr; //!< Increment in the global numbering along the edge + + //! \brief Default constructor. + Edge() { icnod = incr = 0; } + //! \brief Returns \a icnod which then is incremented. + int next(); + }; + +public: + //! \brief Struct with data for definition of global node numbers of a patch. + struct BlockNodes + { + int ibnod[4]; //!< Vertex nodes + Edge edges[4]; //!< Edge nodes + int iinod; //!< Global node number of the first interior node + int inc[2]; //!< Increment in global node numbering in each direction + int nnodI; //!< Number of nodes in parameter direction I + int indxI; //!< Running node index in the local I-direction + + //! \brief Default constructor. + BlockNodes() { iinod = inc[0] = inc[1] = 0; indxI = 1; } + //! \brief Returns \a iinod which then is incremented. + int next(); + }; + +#endif + +public: + + //! \brief Constructor creating an instance by reading the given file. + ASMu2D(const char* fName = 0, unsigned char n_s = 2, unsigned char n_f = 2); + //! \brief Constructor creating an instance by reading the given input stream. + ASMu2D(std::istream& is, unsigned char n_s = 2, unsigned char n_f = 2); + //! \brief Empty destructor. + virtual ~ASMu2D() {} + + + // Methods to access data + // ============================ + LR::LRSplineSurface* getSurface() { return lrspline; } + + // Methods for model generation + // ============================ + + //! \brief Generates the finite element topology data for the patch. + //! \details The data generated are the element-to-node connectivity array, + //! the node-to-IJ-index array, as well as global node and element numbers. + virtual bool generateFEMTopology(); + + //! \brief Clears the contents of the patch, making it empty. + virtual void clear(); + + //! \brief Returns the global coordinates for the given node. + //! \param[in] inod 1-based node index local to current patch + virtual Vec3 getCoord(size_t inod) const; + + //! \brief Creates an instance by reading the given input stream. + bool read(std::istream&); + //! \brief Writes the geometry of the SplineSurface object to given stream. + virtual bool write(std::ostream&, int = 0) const; + +#if 0 + //! \brief Assigns new global node numbers for all nodes of the patch. + //! \param nodes Object with global nodes numbers to assign to this patch + //! \param[in] basis Which basis to assign node numbers for in mixed methods + //! + //! \details The global node numbers generated by \a generateFEMTopology are + //! non-unique in the sense that a node that is shared by two (or more) + //! patches along a common interface has a different number in each patch. + //! This method therefore assigns a new global number to each node in the + //! patch. The data provided through the \a nodes argument is sufficient + //! to determine the unique global number under the assumption that they + //! are ordered in the sequence determined by the local orientation of the + //! patch and its edges. + bool assignNodeNumbers(); +#endif + + //! \brief Progressively refine the diagonal until at least the LR-spline + //! contains at least \a minBasisfunctions basisfunctions + //! \param[in] minBasisfunctions lower bound on number of resulting basis functions + bool diagonalRefine(int minBasisfunctions); + //! \brief Progressively refine the lower left corner until at least the LR-spline + //! contains at least \a minBasisfunctions basisfunctions + //! \param[in] minBasisfunctions lower bound on number of resulting basis functions + bool cornerRefine(int minBasisfunctions); + //! \brief Inserts one (global)line at a time until the LR-spline contains at least + //! \a minBasisfunctions basisfunctions + //! \param[in] minBasisfunctions lower bound on number of resulting basis functions + bool uniformRefine(int minBasisfunctions); + + // Various methods for preprocessing of boundary conditions and patch topology + // =========================================================================== + + //! \brief Constrains all DOFs on a given boundary edge. + //! \param[in] dir Parameter direction defining the edge to constrain + //! \param[in] dof Which DOFs to constrain at each node on the edge + //! \param[in] code Inhomogeneous dirichlet condition code + void constrainEdge(int dir, int dof = 123, int code = 0); + + //! \brief Constrains a corner node identified by the two parameter indices. + //! \param[in] I Parameter index in u-direction + //! \param[in] J Parameter index in v-direction + //! \param[in] dof Which DOFs to constrain at the node + //! \param[in] code Inhomogeneous dirichlet condition code + //! + //! \details The sign of the two indices is used to define whether we want + //! the node at the beginning or the end of that parameter direction. + //! The magnitude of the indices are not used. + void constrainCorner(int I, int J, int dof = 123, int code = 0); + //! \brief Constrains a node identified by two relative parameter values. + //! \param[in] xi Parameter in u-direction + //! \param[in] eta Parameter in v-direction + //! \param[in] dof Which DOFs to constrain at the node + //! \param[in] code Inhomogeneous dirichlet condition code + //! + //! \details The parameter values have to be in the domain [0.0,1.0], where + //! 0.0 means the beginning of the domain and 1.0 means the end. For values + //! in between, the actual index is taken as the integer value closest to + //! \a r*n, where \a r denotes the given relative parameter value, + //! and \a n is the number of nodes along that parameter direction. + void constrainNode(double xi, double eta, int dof = 123, int code = 0); + +// More multipatch stuff +#if 0 + + //! \brief Connects all matching nodes on two adjacent boundary edges. + //! \param[in] edge Local edge index of this patch, in range [1,4] + //! \param neighbor The neighbor patch + //! \param[in] nedge Local edge index of neighbor patch, in range [1,4] + //! \param[in] revers Indicates whether the two edges have opposite directions + virtual bool connectPatch(int edge, ASMu2D& neighbor, int nedge, + bool revers = false); + + //! \brief Makes two opposite boundary edges periodic. + //! \param[in] dir Parameter direction defining the periodic edges + //! \param[in] basis Which basis to connect (mixed methods), 0 means both + //! \param[in] master 1-based index of the first master node in this basis + virtual void closeEdges(int dir, int basis = 0, int master = 1); + +#endif + + // Methods for integration of finite element quantities. + // These are the main computational methods of the ASM class hierarchy. + // ==================================================================== + + //! \brief Evaluates an integral over the interior patch domain. + //! \param integrand Object with problem-specific data and methods + //! \param glbInt The integrated quantity + //! \param[in] time Parameters for nonlinear/time-dependent simulations + //! \param locInt Vector of element-wise contributions to \a glbInt + virtual bool integrate(Integrand& integrand, + GlobalIntegral& glbInt, const TimeDomain& time, + const LintegralVec& locInt = LintegralVec()); + + //! \brief Evaluates a boundary integral over a patch edge. + //! \param integrand Object with problem-specific data and methods + //! \param[in] lIndex Local index of the boundary edge + //! \param glbInt The integrated quantity + //! \param[in] time Parameters for nonlinear/time-dependent simulations + //! \param locInt Vector of element-wise contributions to \a glbInt + virtual bool integrate(Integrand& integrand, int lIndex, + GlobalIntegral& glbInt, const TimeDomain& time, + const LintegralVec& locInt = LintegralVec()); + + + // Post-processing methods + // ======================= + + //! \brief Evaluates the geometry at a specified point. + //! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point + //! \param[out] param The (u,v) parameters of the point in knot-span domain + //! \param[out] X The Cartesian coordinates of the point + //! \return Local node number within the patch that matches the point, if any + //! \return 0 if no node (control point) matches this point + virtual int evalPoint(const double* xi, double* param, Vec3& X) const; + +// we'll figure out the postprocessing later +#if 0 + //! \brief Creates a quad element model of this patch for visualization. + //! \param[out] grid The generated quadrilateral grid + //! \param[in] npe Number of visualization nodes over each knot span + //! \note The number of element nodes must be set in \a grid on input. + virtual bool tesselate(ElementBlock& grid, const int* npe) const; + + //! \brief Evaluates the primary solution field at all visualization points. + //! \param[out] sField Solution field + //! \param[in] locSol Solution vector in DOF-order + //! \param[in] npe Number of visualization nodes over each knot span + virtual bool evalSolution(Matrix& sField, const Vector& locSol, + const int* npe) const; + + //! \brief Evaluates the primary solution field at the given points. + //! \param[out] sField Solution field + //! \param[in] locSol Solution vector local to current patch + //! \param[in] gpar Parameter values of the result sampling points + //! \param[in] regular Flag indicating how the sampling points are defined + //! + //! \details When \a regular is \e true, it is assumed that the parameter + //! value array \a gpar forms a regular tensor-product point grid of dimension + //! \a gpar[0].size() \a X \a gpar[1].size(). + //! Otherwise, we assume that it contains the \a u and \a v parameters + //! directly for each sampling point. + virtual bool evalSolution(Matrix& sField, const Vector& locSol, + const RealArray* gpar, bool regular = true) const; + + //! \brief Evaluates the secondary solution field at all visualization points. + //! \param[out] sField Solution field + //! \param[in] integrand Object with problem-specific data and methods + //! \param[in] npe Number of visualization nodes over each knot span + //! \param[in] project Flag indicating if the projected solution is wanted + //! + //! \details The secondary solution is derived from the primary solution, + //! which is assumed to be stored within the \a integrand for current patch. + //! If \a npe is NULL, the solution is evaluated at the Greville points and + //! then projected onto the spline basis to obtain the control point values, + //! which then are returned through \a sField. + //! If \a npe is not NULL and \a project is \e true, the solution is also + //! projected onto the spline basis, and then evaluated at the \a npe points + virtual bool evalSolution(Matrix& sField, const Integrand& integrand, + const int* npe = 0, bool project = false) const; + + //! \brief Projects the secondary solution field onto the primary basis. + //! \param[in] integrand Object with problem-specific data and methods + Go::SplineSurface* projectSolution(const Integrand& integrand) const; + //! \brief Projects the secondary solution field onto the primary basis. + //! \param[in] integrand Object with problem-specific data and methods + virtual Go::GeomObject* evalSolution(const Integrand& integrand) const; + + //! \brief Evaluates the secondary solution field at the given points. + //! \param[out] sField Solution field + //! \param[in] integrand Object with problem-specific data and methods + //! \param[in] gpar Parameter values of the result sampling points + //! \param[in] regular Flag indicating how the sampling points are defined + //! + //! \details The secondary solution is derived from the primary solution, + //! which is assumed to be stored within the \a integrand for current patch. + //! When \a regular is \e true, it is assumed that the parameter value array + //! \a gpar forms a regular tensor-product point grid of dimension + //! \a gpar[0].size() \a X \a gpar[1].size(). + //! Otherwise, we assume that it contains the \a u and \a v parameters + //! directly for each sampling point. + virtual bool evalSolution(Matrix& sField, const Integrand& integrand, + const RealArray* gpar, bool regular = true) const; + + //! \brief Calculates parameter values for visualization nodal points. + //! \param[out] prm Parameter values in given direction for all points + //! \param[in] dir Parameter direction (0,1) + //! \param[in] nSegSpan Number of visualization segments over each knot-span + virtual bool getGridParameters(RealArray& prm, int dir, int nSegSpan) const; +#endif + +protected: + + // Internal utility methods + // ======================== + + //! \brief Connects all matching nodes on two adjacent boundary edges. + //! \param[in] edge Local edge index of this patch, in range [1,4] + //! \param neighbor The neighbor patch + //! \param[in] nedge Local edge index of neighbor patch, in range [1,4] + //! \param[in] revers Indicates whether the two edges have opposite directions + //! \param[in] basis Which basis to connect the nodes for (mixed methods) + //! \param[in] slave 0-based index of the first slave node in this basis + //! \param[in] master 0-based index of the first master node in this basis + bool connectBasis(int edge, ASMu2D& neighbor, int nedge, bool revers, + int basis = 1, int slave = 0, int master = 0); + + //! \brief Extracts parameter values of the Gauss points in one direction. + //! \param[out] uGP Parameter values in given direction for all points + //! \param[in] dir Parameter direction (0,1) + //! \param[in] nGauss Number of Gauss points along a knot-span + //! \param[in] iel Element index + //! \param[in] xi Dimensionless Gauss point coordinates [-1,1] + void getGaussPointParameters(Vector& uGP, int dir, int nGauss, + int iel, const double* xi) const; + + //! \brief Calculates parameter values for the Greville points. + //! \param[out] prm Parameter values in given direction for all points + //! \param[in] dir Parameter direction (0,1) + bool getGrevilleParameters(RealArray& prm, int dir) const; + + //! \brief Returns the area in the parameter space for an element. + //! \param[in] iel Element index + 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 + double getParametricLength(int iel, int dir) const; + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes + //! in one element + virtual bool getElementCoordinates(Matrix& X, int iel) const; + //! \brief Returns a matrix with all nodal coordinates within the patch. + //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes + //! in the patch + virtual void getNodalCoordinates(Matrix& X) const; + + //! \brief Establishes matrices with basis functions and 1st derivatives. + static void extractBasis(const Go::BasisDerivsSf& spline, + Vector& N, Matrix& dNdu); + //! \brief Establishes matrices with basis functions, 1st and 2nd derivatives. + static void extractBasis(const Go::BasisDerivsSf2& spline, + Vector& N, Matrix& dNdu, Matrix3D& d2Ndu2); + + //! \brief Auxilliary function for computation of basis function indices. + static void scatterInd(int n1, int n2, int p1, int p2, + const int* start, IntVec& index); + +private: + + +protected: + LR::LRSplineSurface* lrspline; //!< Pointer to the actual spline surface object +}; + +#endif diff --git a/src/ASM/LR/ASMunstruct.C b/src/ASM/LR/ASMunstruct.C new file mode 100644 index 00000000..18e4a4b5 --- /dev/null +++ b/src/ASM/LR/ASMunstruct.C @@ -0,0 +1,34 @@ +// $Id$ +//============================================================================== +//! +//! \file ASMunstruct.C +//! +//! \date December 2010 +//! +//! \author Kjetil Andre Johannessen / SINTEF +//! +//! \brief Base class for unstructured spline-based FE assembly drivers. +//! +//============================================================================== + +#include "ASMunstruct.h" +#include "GoTools/geometry/GeomObject.h" +#include "LRSpline/LRSplineSurface.h" + + +int ASMunstruct::gEl = 0; +int ASMunstruct::gNod = 0; + + +ASMunstruct::ASMunstruct (unsigned char n_p, unsigned char n_s, unsigned char n_f) + : ASMbase(n_p,n_s,n_f) +{ + geo = 0; +} + + +ASMunstruct::~ASMunstruct () +{ + if (geo) delete geo; +} + diff --git a/src/ASM/LR/ASMunstruct.h b/src/ASM/LR/ASMunstruct.h new file mode 100644 index 00000000..b5405048 --- /dev/null +++ b/src/ASM/LR/ASMunstruct.h @@ -0,0 +1,65 @@ +// $Id$ +//============================================================================== +//! +//! \file ASMunstruct.h +//! +//! \date September 2011 +//! +//! \author Kjetil Andre Johannessen / SINTEF +//! +//! \brief Base class for unstructured spline-based FE assembly drivers. +//! +//============================================================================== + +#ifndef _ASM_UNSTRUCT_H +#define _ASM_UNSTRUCT_H + +#include "ASMbase.h" + +namespace Go { + class GeomObject; + class BoundingBox; +} +namespace LR { + class LRSplineSurface; +} + +/*! + \brief Base class for structured spline-based FE assembly drivers. + \details This class contains methods common for structured spline patches. +*/ + +class ASMunstruct : public ASMbase +{ +protected: + //! \brief The constructor sets the space dimensions. + //! \param[in] n_p Number of parameter dimensions + //! \param[in] n_s Number of spatial dimensions + //! \param[in] n_f Number of primary solution fields + ASMunstruct(unsigned char n_p, unsigned char n_s, unsigned char n_f); + +public: + //! \brief The destructor frees the dynamically allocated spline object. + virtual ~ASMunstruct(); + + //! \brief Checks if the patch is empty. + virtual bool empty() const { return geo == 0; } + + //! \brief Defines the numerical integration scheme to use. + //! \param[in] ng Number of Gauss points in each parameter direction + void setGauss(int ng) { nGauss = ng; } + + //! \brief Resets global element and node counters + static void resetNumbering() { gEl = gNod = 0; } + +protected: + LR::LRSplineSurface* geo; //!< Pointer to the actual spline geometry object + + int nGauss; //!< Numerical integration scheme + static int gEl; //!< Global element counter + static int gNod; //!< Global node counter + +}; + +#endif + diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C index 5d8abbc0..e8945546 100644 --- a/src/SIM/SIM2D.C +++ b/src/SIM/SIM2D.C @@ -11,6 +11,9 @@ //! //============================================================================== +#ifdef HAS_LRSPLINE +#include "LR/ASMu2D.h" +#endif #include "SIM2D.h" #include "ASMs2Dmx.h" #include "ASMs2DmxLag.h" @@ -51,6 +54,11 @@ bool SIM2D::parse (char* keyWord, std::istream& is) case Spectral: pch = new ASMs2DSpec(cline,2,nf[0]); break; + #ifdef HAS_LRSPLINE + case LRSpline: + pch = new ASMu2D(cline,2,nf[0]); + break; + #endif default: if (nf[1] > 0) pch = new ASMs2Dmx(cline,2,nf[0],nf[1]); @@ -120,19 +128,21 @@ bool SIM2D::parse (char* keyWord, std::istream& is) if (!oneBasedIdx) { - // We always require the node numbers to be 1-based - for (i = 0; i < 4; i++) ++n.ibnod[i]; - for (i = 0; i < 4; i++) ++n.edges[i].icnod; - ++n.iinod; + // We always require the node numbers to be 1-based + for (i = 0; i < 4; i++) ++n.ibnod[i]; + for (i = 0; i < 4; i++) ++n.edges[i].icnod; + ++n.iinod; } if (isn.good() && pid > 0) - if (!static_cast(myModel[pid-1])->assignNodeNumbers(n)) - { - std::cerr <<" *** SIM2D::parse: Failed to assign node numbers" - <<" for patch "<< patch << std::endl; - return false; - } + { + if (!static_cast(myModel[pid-1])->assignNodeNumbers(n)) + { + std::cerr <<" *** SIM2D::parse: Failed to assign node numbers" + <<" for patch "<< patch << std::endl; + return false; + } + } } } @@ -431,40 +441,79 @@ bool SIM2D::addConstraint (int patch, int lndx, int ldim, int dirs, int code) return constrError("patch index ",patch); std::cout <<"\tConstraining P"<< patch - << (ldim == 0 ? " V" : " E") << lndx <<" in direction(s) "<< dirs; + << (ldim == 0 ? " V" : " E") << lndx <<" in direction(s) "<< dirs; if (code) std::cout <<" code = "<< code <<" "; - ASMs2D* pch = static_cast(myModel[patch-1]); - switch (ldim) - { - case 0: // Vertex constraints - switch (lndx) - { - case 1: pch->constrainCorner(-1,-1,dirs,code); break; - case 2: pch->constrainCorner( 1,-1,dirs,code); break; - case 3: pch->constrainCorner(-1, 1,dirs,code); break; - case 4: pch->constrainCorner( 1, 1,dirs,code); break; - default: std::cout << std::endl; - return constrError("vertex index ",lndx); - } - break; - - case 1: // Edge constraints - switch (lndx) - { - case 1: pch->constrainEdge(-1,dirs,code); break; - case 2: pch->constrainEdge( 1,dirs,code); break; - case 3: pch->constrainEdge(-2,dirs,code); break; - case 4: pch->constrainEdge( 2,dirs,code); break; - default: std::cout << std::endl; - return constrError("edge index ",lndx); - } - break; - - default: - std::cout << std::endl; - return constrError("local dimension switch ",ldim); - } + if(discretization == LRSpline ) { + #ifdef HAS_LRSPLINE + ASMu2D* pch = static_cast(myModel[patch-1]); + switch (ldim) + { + case 0: // Vertex constraints + switch (lndx) + { + case 1: pch->constrainCorner(-1,-1,dirs,code); break; + case 2: pch->constrainCorner( 1,-1,dirs,code); break; + case 3: pch->constrainCorner(-1, 1,dirs,code); break; + case 4: pch->constrainCorner( 1, 1,dirs,code); break; + default: std::cout << std::endl; + return constrError("vertex index ",lndx); + } + break; + + case 1: // Edge constraints + switch (lndx) + { + case 1: pch->constrainEdge(-1,dirs,code); break; + case 2: pch->constrainEdge( 1,dirs,code); break; + case 3: pch->constrainEdge(-2,dirs,code); break; + case 4: pch->constrainEdge( 2,dirs,code); break; + default: std::cout << std::endl; + return constrError("edge index ",lndx); + } + break; + + default: + std::cout << std::endl; + return constrError("local dimension switch ",ldim); + } + #else + std::cerr << "Error: No LR-spline library detected\n"; + return constrError("discretization: LR-spline ", discretization); + #endif + } else { + ASMs2D* pch = static_cast(myModel[patch-1]); + switch (ldim) + { + case 0: // Vertex constraints + switch (lndx) + { + case 1: pch->constrainCorner(-1,-1,dirs,code); break; + case 2: pch->constrainCorner( 1,-1,dirs,code); break; + case 3: pch->constrainCorner(-1, 1,dirs,code); break; + case 4: pch->constrainCorner( 1, 1,dirs,code); break; + default: std::cout << std::endl; + return constrError("vertex index ",lndx); + } + break; + + case 1: // Edge constraints + switch (lndx) + { + case 1: pch->constrainEdge(-1,dirs,code); break; + case 2: pch->constrainEdge( 1,dirs,code); break; + case 3: pch->constrainEdge(-2,dirs,code); break; + case 4: pch->constrainEdge( 2,dirs,code); break; + default: std::cout << std::endl; + return constrError("edge index ",lndx); + } + break; + + default: + std::cout << std::endl; + return constrError("local dimension switch ",ldim); + } + } return true; } @@ -494,6 +543,11 @@ void SIM2D::readPatches (std::istream& isp) case Spectral: pch = new ASMs2DSpec(isp,2,nf[0]); break; + case LRSpline: + #ifdef HAS_LRSPLINE + pch = new ASMu2D(isp,2,nf[0]); + #endif + break; default: if (nf[1] > 0) pch = new ASMs2Dmx(isp,2,nf[0],nf[1]); diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index fa3511c0..01e0454f 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -492,7 +492,7 @@ protected: public: //! \brief Enum defining the available discretization methods. - enum Discretization { Spline, Lagrange, Spectral }; + enum Discretization { Spline, Lagrange, Spectral, LRSpline }; static Discretization discretization; //!< Spatial discretization option