From b335397a03759aa5feaa12cdb39debfbe01b0049 Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Tue, 22 Dec 2020 11:06:13 +0100 Subject: [PATCH 01/11] Missing decorators to get it compiling with lrspline-python and IFEM --- include/LRSpline/Basisfunction.h | 2 +- src/Basisfunction.cpp | 2 +- src/MeshRectangle.cpp | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/include/LRSpline/Basisfunction.h b/include/LRSpline/Basisfunction.h index 75e03f9..2e42097 100644 --- a/include/LRSpline/Basisfunction.h +++ b/include/LRSpline/Basisfunction.h @@ -170,7 +170,7 @@ class Basisfunction : public Streamable { std::vector::const_iterator cp() const { return controlpoint_.begin(); }; double cp(int i) const { return controlpoint_[i]; }; double w() const { return weight_; }; - double integral(Element *el) const ; + double integral(const Element *el) const ; long hashCode() const ; diff --git a/src/Basisfunction.cpp b/src/Basisfunction.cpp index 10d1956..66f89d5 100644 --- a/src/Basisfunction.cpp +++ b/src/Basisfunction.cpp @@ -848,7 +848,7 @@ void Basisfunction::normalize(int pardir, double parmin, double parmax) { /************************************************************************************************************************//** * \brief computes the integral of this basisfunction over a particular element. ***************************************************************************************************************************/ -double Basisfunction::integral(Element* el) const { +double Basisfunction::integral(const Element* el) const { /* This function is built on two observations: * 1. Johannessen K.A (https://doi.org/10.1016/j.cma.2016.04.030), Equation (5) section 2.1 (integration of splines) * \int N_i,p,t = (t_{i+p+1}-t_i)/(p+1) \sum_j N_j,p+1,T, where T=(t0,t1,t2,..,tn,tn) diff --git a/src/MeshRectangle.cpp b/src/MeshRectangle.cpp index 52223f8..04dafa6 100644 --- a/src/MeshRectangle.cpp +++ b/src/MeshRectangle.cpp @@ -3,6 +3,7 @@ #include "LRSpline/Element.h" #include "LRSpline/Basisfunction.h" #include +#include namespace LR { From bb75482a14d26a3b57f140720a7b7d44e60c26e9 Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Tue, 22 Dec 2020 12:12:32 +0100 Subject: [PATCH 02/11] Added: counting number of refinements per element --- Apps/TestLevelCounting.cpp | 141 ++++++++++++++++++++++++++++ CMakeLists.txt | 8 ++ include/LRSpline/Element.h | 19 +++- src/Element.cpp | 33 ++++++- test/TestLevelCounting/testSurf.reg | 9 ++ test/TestLevelCounting/testVol.reg | 8 ++ 6 files changed, 208 insertions(+), 10 deletions(-) create mode 100644 Apps/TestLevelCounting.cpp create mode 100644 test/TestLevelCounting/testSurf.reg create mode 100644 test/TestLevelCounting/testVol.reg diff --git a/Apps/TestLevelCounting.cpp b/Apps/TestLevelCounting.cpp new file mode 100644 index 0000000..9ce2d5a --- /dev/null +++ b/Apps/TestLevelCounting.cpp @@ -0,0 +1,141 @@ +#include +#include +#include +#include +#include +#include +#include "LRSpline/LRSplineSurface.h" +#include "LRSpline/LRSplineVolume.h" +#include "LRSpline/MeshRectangle.h" +#include "LRSpline/Meshline.h" +#include "LRSpline/Profiler.h" +#include "LRSpline/Element.h" + +using namespace Go; +using namespace LR; +using namespace std; + +int main(int argc, char **argv) { + + // set default parameter values + int p1 = 3; + int p2 = 3; + int p3 = 3; + int n1 = 9; + int n2 = 9; + int n3 = 9; + int dim = 4; + bool rat = false; + bool vol = false; + string parameters(" parameters: \n" \ + " -p1 polynomial ORDER (degree+1) in first parametric direction\n" \ + " -p2 polynomial order in second parametric direction\n" \ + " -p3 polynomial order in third parametric direction\n" \ + " -n1 number of basis functions in first parametric direction\n" \ + " -n2 number of basis functions in second parametric direction\n" \ + " -n3 number of basis functions in third parametric direction\n" \ + " -dim dimension of the controlpoints\n" \ + " -vol enforce trivariate volumetric test case\n" \ + " -help display (this) help screen\n"); + + // read input + for(int i=1; i knot_u(n1 + p1); + std::vector knot_v(n2 + p2); + std::vector knot_w(n3 + p3); + for(int i=0; in1) ? n1-p1+1 : i-p1+1; + for(int i=0; in2) ? n2-p2+1 : i-p2+1; + for(int i=0; in3) ? n3-p3+1 : i-p3+1; + + // create a list of random control points (all between 0.1 and 1.1) + int nCP = (vol) ? n1*n2*n3 : n1*n2; + nCP *= (dim+rat); + std::vector cp(nCP); + int k=0; + for(int i=0; i0 + + if(vol) { + lr = lrv = new LRSplineVolume(n1, n2, n3, p1, p2, p3, knot_u.begin(), knot_v.begin(), knot_w.begin(), cp.begin(), dim, rat); + } else { + lr = lrs = new LRSplineSurface(n1, n2, p1, p2, knot_u.begin(), knot_v.begin(), cp.begin(), dim, rat); + } + // refine the lower-left corner 4 times + for(int i=0; i<4; i++) { + lr->generateIDs(); + vector corner(lr->nVariate(), 0.00001); + Element* el = lr->getElement(lr->getElementContaining(corner)); + vector refI(0); + for(auto b : el->support()) refI.push_back(b->getId()); + lr->refineBasisFunction(refI); + } + lr->generateIDs(); + + bool allOK = true; + for(auto el : lr->getAllElements()) { + double du = el->umax() - el->umin(); + int expectedLevel = round(log2(1/du)); + cout << "Element #" << el->getId() << ": Length=" << du << " level " << el->getLevel(0); + for(int i=1; inVariate(); i++) cout << ", " << el->getLevel(i); + cout << " (expects lvl " << expectedLevel << ")"; + bool ok = true; + for(int i=0; inVariate(); i++) + ok &= expectedLevel == el->getLevel(i); + if(ok) cout << " - OK" << endl; + else cout << " - ASSERTION FAIL" << endl; + allOK &= ok; + } + cout << "====================================================" << endl; + if(allOK) cout << "All assertions passed" << endl; + else cout << "FAILED TEST" << endl; + cout << "----------------------------------------------------" << endl; +} + diff --git a/CMakeLists.txt b/CMakeLists.txt index 7a6ba42..22dd8c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -172,6 +172,9 @@ ENDIF(HAS_GOTOOLS) ADD_EXECUTABLE(TopologyRefinement ${PROJECT_SOURCE_DIR}/Apps/TopologyRefinement.cpp) TARGET_LINK_LIBRARIES(TopologyRefinement LRSpline ${DEPLIBS}) +ADD_EXECUTABLE(TestLevelCounting ${PROJECT_SOURCE_DIR}/Apps/TestLevelCounting.cpp) +TARGET_LINK_LIBRARIES(TestLevelCounting LRSpline ${DEPLIBS}) + IF(HAS_BOOST) ADD_EXECUTABLE(LinearIndep ${PROJECT_SOURCE_DIR}/Apps/LinearIndep.cpp) TARGET_LINK_LIBRARIES(LinearIndep LRSpline ${DEPLIBS}) @@ -231,6 +234,11 @@ FOREACH(TESTFILE ${REGRESSION_TESTFILES}) ADD_TEST(${TESTFILE} ${PROJECT_SOURCE_DIR}/test/regtest.sh "${CMAKE_BINARY_DIR}/${EXECUTABLE_OUTPUT_PATH}/TopologyRefinement" "${TESTFILE}") ENDFOREACH() +FILE(GLOB_RECURSE REGRESSION_TESTFILES "${PROJECT_SOURCE_DIR}/test/TestLevelCounting/*.reg") +FOREACH(TESTFILE ${REGRESSION_TESTFILES}) + ADD_TEST(${TESTFILE} ${PROJECT_SOURCE_DIR}/test/regtest.sh "${CMAKE_BINARY_DIR}/${EXECUTABLE_OUTPUT_PATH}/TestLevelCounting" "${TESTFILE}") +ENDFOREACH() + FILE(GLOB_RECURSE REGRESSION_TESTFILES "${PROJECT_SOURCE_DIR}/test/Integrals/*.reg") FOREACH(TESTFILE ${REGRESSION_TESTFILES}) ADD_TEST(${TESTFILE} ${PROJECT_SOURCE_DIR}/test/regtest.sh "${CMAKE_BINARY_DIR}/${EXECUTABLE_OUTPUT_PATH}/testIntegral" "${TESTFILE}") diff --git a/include/LRSpline/Element.h b/include/LRSpline/Element.h index 69f547f..d35689b 100644 --- a/include/LRSpline/Element.h +++ b/include/LRSpline/Element.h @@ -36,7 +36,8 @@ class Element : public Streamable { std::copy(lowerLeft, lowerLeft + dim, min.begin()); std::copy(upperRight, upperRight + dim, max.begin()); id_ = -1; - overloadCount = 0; + overloadCount_= 0; + level_ = std::vector(dim, 0); } Element(std::vector &lowerLeft, std::vector &upperRight); void removeSupportFunction(Basisfunction *f); @@ -79,9 +80,14 @@ class Element : public Streamable { void setVmax(double v) { max[1] = v; }; bool isOverloaded() const; - void resetOverloadCount() { overloadCount = 0; } - int incrementOverloadCount() { return overloadCount++; } - int getOverloadCount() const { return overloadCount; } + void resetOverloadCount() { overloadCount_ = 0; } + int incrementOverloadCount() { return overloadCount_++; } + int getOverloadCount() const { return overloadCount_; } + + //! set the level (number of refinements) of this element (counting different in all directions) + void setLevel(int direction, int new_lvl) { if(direction < min.size()) level_[direction] = new_lvl; } + //! get the level (number of refinements) of this element (counting different in all directions) + int getLevel(int direction) const { return level_[direction]; } void updateBasisPointers(std::vector &basis) ; @@ -96,7 +102,10 @@ class Element : public Streamable { HashSet support_; std::vector support_ids_; // temporary storage for the read() method only - int overloadCount ; + int overloadCount_ ; + // number of refinements for this particular element (i.e. hierarchical B-spline refinement) + // This is counting separately in each parametric direction + std::vector level_ ; }; diff --git a/src/Element.cpp b/src/Element.cpp index c23a516..1d0152f 100644 --- a/src/Element.cpp +++ b/src/Element.cpp @@ -12,8 +12,9 @@ namespace LR { * \brief Default constructor ***************************************************************************************************************************/ Element::Element() { - id_ = -1; - overloadCount = 0; + id_ = -1; + overloadCount_ = 0; + level_ = std::vector(0); } /************************************************************************************************************************//** @@ -22,7 +23,8 @@ Element::Element() { ***************************************************************************************************************************/ Element::Element(int dim) { id_ = -1; - overloadCount = 0; + overloadCount_= 0; + level_ = std::vector(dim,0); min.resize(dim); max.resize(dim); } @@ -43,7 +45,8 @@ Element::Element(double start_u, double start_v, double stop_u, double stop_v) { max[0] = stop_u ; max[1] = stop_v ; id_ = -1; - overloadCount = 0; + overloadCount_ = 0; + level_ = std::vector(2,0); } /************************************************************************************************************************//** @@ -71,6 +74,7 @@ Element* Element::copy() { returnvalue->id_ = this->id_; returnvalue->min = this->min; // it seems that the default vector operator= thing takes a deep copy returnvalue->max = this->max; + returnvalue->level_ = this->level_; for(Basisfunction* b : support_) returnvalue->support_ids_.push_back(b->getId()); @@ -97,6 +101,9 @@ Element* Element::split(int splitDim, double par_value) { max[splitDim] = par_value; // old element should stop at par_value newElement = new Element(min.size(), newMin.begin(), newMax.begin()); + level_[splitDim]++; + for(uint dim=0; dimsetLevel(dim, level_[dim]); for(Basisfunction *b : support_) if(b->addSupport(newElement)) // tests for overlapping as well @@ -156,6 +163,7 @@ void Element::read(std::istream &is) { ASSERT_NEXT_CHAR(':'); min.resize(dim); max.resize(dim); + level_.resize(dim); ASSERT_NEXT_CHAR('('); is >> min[0]; @@ -187,6 +195,19 @@ void Element::read(std::istream &is) { nextChar = is.peek(); } ASSERT_NEXT_CHAR('}'); + nextChar = is.peek(); + if(nextChar == 'l') { // after v.1.6.0 (level specified) + ASSERT_NEXT_CHAR('l'); + ASSERT_NEXT_CHAR('v'); + ASSERT_NEXT_CHAR('l'); + is >> level_[0]; + for(int i=1; i> level_[i]; + } + } else { // else level undefined, default to all 0's + level_ = std::vector(dim,0); + } } #undef ASSERT_NEXT_CHAR @@ -210,7 +231,9 @@ void Element::write(std::ostream &os) const { os << b->getId(); isFirst = false; } - os << "}"; + os << "} lvl " << level_[0]; + for(uint i=1; i Date: Tue, 22 Dec 2020 13:43:05 +0100 Subject: [PATCH 03/11] Added force max T-joints for volumetric mesh --- include/LRSpline/Element.h | 4 ++++ include/LRSpline/LRSplineVolume.h | 4 ++-- src/Element.cpp | 39 ++++++++++++++++++++++++++++++ src/LRSplineVolume.cpp | 40 +++++++++++++++++++++++++++++++ 4 files changed, 85 insertions(+), 2 deletions(-) diff --git a/include/LRSpline/Element.h b/include/LRSpline/Element.h index d35689b..8f070da 100644 --- a/include/LRSpline/Element.h +++ b/include/LRSpline/Element.h @@ -89,6 +89,10 @@ class Element : public Streamable { //! get the level (number of refinements) of this element (counting different in all directions) int getLevel(int direction) const { return level_[direction]; } + // topological operators + bool touches(const Element* el) const ; + std::vector neighbours() ; + void updateBasisPointers(std::vector &basis) ; virtual void read(std::istream &is); diff --git a/include/LRSpline/LRSplineVolume.h b/include/LRSpline/LRSplineVolume.h index 596c6bd..80f66de 100644 --- a/include/LRSpline/LRSplineVolume.h +++ b/include/LRSpline/LRSplineVolume.h @@ -91,11 +91,11 @@ class LRSplineVolume : public LRSpline { /* MeshRectangle* insert_const_u_edge(double u, double start_v, double stop_v, int multiplicity=1); MeshRectangle* insert_const_v_edge(double v, double start_u, double stop_u, int multiplicity=1); - void aPosterioriFixes() ; void closeGaps( std::vector* newLines=NULL); - void enforceMaxTjoints( std::vector* newLines=NULL); void enforceMaxAspectRatio(std::vector* newLines=NULL); */ + void aPosterioriFixes() ; + void enforceMaxTjoints(); bool enforceIsotropic(); // linear independence methods diff --git a/src/Element.cpp b/src/Element.cpp index 1d0152f..d4dde5e 100644 --- a/src/Element.cpp +++ b/src/Element.cpp @@ -3,6 +3,7 @@ #include "LRSpline/Meshline.h" #include "LRSpline/Basisfunction.h" #include +#include typedef unsigned int uint; @@ -124,6 +125,44 @@ Element* Element::split(int splitDim, double par_value) { return newElement; } +/************************************************************************************************************************//** + * \brief returns true if the element shares a side with the current element (corners does NOT count, lines in 3D does NOT count) + ***************************************************************************************************************************/ +bool Element::touches(const Element* el) const { + for(uint dim=0; dimgetParmin(dim) == this->getParmax(dim) || el->getParmax(dim) == this->getParmin(dim)) { + bool overlaps = true; + for(uint i=0; igetParmax(i) > this->getParmin(i)); + overlaps = overlaps && (el->getParmin(i) < this->getParmax(i)); + } + if(overlaps) return true; + } + } + return false; +} + +/************************************************************************************************************************//** + * \brief returns all Elements which share a side (line in 2D, surface in 3D) with this element + ***************************************************************************************************************************/ +std::vector Element::neighbours() { + std::vector result; + + // look at the (element-support) of all functions living on this element + // This is the "extended support" of this element, and the true element + // neigbours are part of this set + std::set potentialNeighbours; + for(Basisfunction* b : this->support()) + for(Element* e : b->support()) + potentialNeighbours.insert(e); + for(Element* el : potentialNeighbours) + if(el->touches(this)) + result.push_back(el); + + return result; +} + /************************************************************************************************************************//** * \brief returns the parametric midpoint of the element, i.e. [ (umin()+umax())/2, (vmin()+vmax())/2] for 2D elements ***************************************************************************************************************************/ diff --git a/src/LRSplineVolume.cpp b/src/LRSplineVolume.cpp index b78d21d..652e1f7 100644 --- a/src/LRSplineVolume.cpp +++ b/src/LRSplineVolume.cpp @@ -1122,6 +1122,46 @@ bool LRSplineVolume::matchParametricEdge(parameterEdge edge, LRSplineVolume *oth } +void LRSplineVolume::aPosterioriFixes() { + int nFunc; + do { + nFunc = basis_.size(); +// if(doCloseGaps_) +// this->closeGaps(newLines); + if(maxTjoints_ > 0) + this->enforceMaxTjoints(); +// if(doAspectRatioFix_) +// this->enforceMaxAspectRatio(newLines); + } while(nFunc != basis_.size()); +} + + +void LRSplineVolume::enforceMaxTjoints() { + bool someFix = true; + while(someFix) { + std::vector newRefinement; + someFix = false; + for(Element* el : element_) { + for(Element* neigh : el->neighbours()) { + if(neigh->getLevel(0) > el->getLevel(0)+1) { + for(int i=0; i<3; i++) { + double start[3] = {el->umin(), el->vmin(), el->wmin()}; + double stop[3] = {el->umax(), el->vmax(), el->wmax()}; + start[i] = (el->getParmin(i)+el->getParmax(i)) / 2.0; + stop[i] = (el->getParmin(i)+el->getParmax(i)) / 2.0; + newRefinement.push_back(new MeshRectangle(start,stop, refKnotlineMult_)); + someFix = true; + } + break; + } + } + } + for(auto m : newRefinement) this->insert_line(m); + } +} + + + /************************************************************************************************************************//** * \brief Enforces all elements to be of equal size in all 3 directions * \return true if some elements were bisected From ac4d41affa2a2721db99546faca78d7d54605763 Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Tue, 22 Dec 2020 13:43:55 +0100 Subject: [PATCH 04/11] bump version 1.6.1 --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 22dd8c7..a846dcb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ PROJECT(LRSpline) CMAKE_MINIMUM_REQUIRED(VERSION 2.6) SET(LRSpline_VERSION_MAJOR 1) SET(LRSpline_VERSION_MINOR 6) -SET(LRSpline_VERSION_PATCH 0) +SET(LRSpline_VERSION_PATCH 1) SET(LRSpline_VERSION ${LRSpline_VERSION_MAJOR}.${LRSpline_VERSION_MINOR}.${LRSpline_VERSION_PATCH}) IF(NOT WIN32) SET(TIME_LRSPLINE 1) From e7a1954ab58b30d0c51b48b53173dd9a25617231 Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Tue, 22 Dec 2020 14:28:11 +0100 Subject: [PATCH 05/11] Updated LinearIndep test to include volumetric inputs --- Apps/LinearIndep.cpp | 56 +++++++++++++++++++++++++++--------------- src/LRSplineVolume.cpp | 23 +++++++++++++++++ 2 files changed, 59 insertions(+), 20 deletions(-) diff --git a/Apps/LinearIndep.cpp b/Apps/LinearIndep.cpp index 59fc87b..3271b2b 100644 --- a/Apps/LinearIndep.cpp +++ b/Apps/LinearIndep.cpp @@ -5,6 +5,7 @@ #include #include #include "LRSpline/LRSplineSurface.h" +#include "LRSpline/LRSplineVolume.h" #include "LRSpline/Profiler.h" #include "LRSpline/Element.h" #include "LRSpline/Meshline.h" @@ -92,7 +93,9 @@ int main(int argc, char **argv) { exit(3); } - LRSplineSurface *lr; + LRSpline *lr; + LRSplineSurface *lrs; + LRSplineVolume *lrv; ifstream inputfile; inputfile.open(inputFileName); @@ -100,11 +103,20 @@ int main(int argc, char **argv) { cerr << "Error: could not open file " << inputFileName << endl; exit(3); } - lr = new LRSplineSurface(); - inputfile >> *lr; + + char buffer[512]; + inputfile.getline(buffer, 512); // peek the first line to figure out if it's an LRSpline or a GoTools spline + inputfile.seekg(ios_base::beg); + if(strncmp(buffer, "# LRSPLINE VOLUME",17)==0) { + lr = lrv = new LRSplineVolume(); + inputfile >> *lrv; + } else { + lr = lrs = new LRSplineSurface(); + inputfile >> *lrs; + } if(isInteger) - lr->makeIntegerKnots(); + lrs->makeIntegerKnots(); bool isLinearIndep = true; if(refineFileName != NULL) { @@ -132,20 +144,20 @@ int main(int argc, char **argv) { /* setting up for refinement */ if(maxTjoints > 0) - lr->setMaxTjoints(maxTjoints); + lrs->setMaxTjoints(maxTjoints); if(maxAspectRatio > 0) - lr->setMaxAspectRatio(maxAspectRatio); - lr->setCloseGaps(closeGaps); + lrs->setMaxAspectRatio(maxAspectRatio); + lrs->setCloseGaps(closeGaps); cout << "calling LRSplineSurface::refine(...)\n"; LRSplineSurface *lr_original = NULL; if(one_by_one) - lr_original = lr->copy(); + lr_original = lrs->copy(); cout << setprecision(16); vector *newLines = new vector(); - // lr->refine(sorted_list, beta, mult, strat, symmetry, newLines); + // lrs->refine(sorted_list, beta, mult, strat, symmetry, newLines); cout << "Number of new mesh lines: " << newLines->size() << endl; for(unsigned int i=0; isize(); i++) { newLines->at(i)->writeMore(cout); @@ -205,7 +217,7 @@ int main(int argc, char **argv) { #endif if( ! isLinearIndep) { printf("Nelements = %5d Nbasis = %5d \n", lr_original->nElements(), lr_original->nBasisFunctions()); - lr = lr_original; + lrs = lr_original; isLinearIndep = false; break; } @@ -214,40 +226,44 @@ int main(int argc, char **argv) { } if(!one_by_one) { - if(floatingPointCheck) - isLinearIndep = lr->isLinearIndepByFloatingPointMappingMatrix(verbose); - else if(overload) + if(floatingPointCheck) { + cout << "Testing for independence by verifying full-rank of mapping matrix to tensor mesh (using floating point numbers)" << endl; + isLinearIndep = lrs->isLinearIndepByFloatingPointMappingMatrix(verbose); + } else if(overload) { + cout << "Testing for independence by overloaded elements" << endl; isLinearIndep = lr->isLinearIndepByOverloading(verbose); #ifdef HAS_BOOST - else + } else { + cout << "Testing for independence by verifying full-rank of mapping matrix to tensor mesh (exact edition)" << endl; isLinearIndep = lr->isLinearIndepByMappingMatrix(verbose); #endif + } } if(dumpFile) { ofstream meshfile; meshfile.open("mesh.eps"); - lr->writePostscriptMesh(meshfile); + lrs->writePostscriptMesh(meshfile); meshfile.close(); ofstream functionfile; functionfile.open("functions.eps"); - lr->writePostscriptFunctionSpace(functionfile); + lrs->writePostscriptFunctionSpace(functionfile); functionfile.close(); ofstream domainfile; domainfile.open("domain.eps"); - lr->writePostscriptElements(domainfile, 10, 10); + lrs->writePostscriptElements(domainfile, 10, 10); domainfile.close(); ofstream controlmesh; controlmesh.open("controlmesh.eps"); - lr->writePostscriptMeshWithControlPoints(controlmesh, 10, 10); + lrs->writePostscriptMeshWithControlPoints(controlmesh, 10, 10); controlmesh.close(); ofstream lrfile; lrfile.open("lrspline.txt"); - lrfile << *lr << endl; + lrfile << *lrs << endl; lrfile.close(); cout << endl; @@ -262,7 +278,7 @@ int main(int argc, char **argv) { #ifdef HAS_BOOST if(dumpNullSpace) { vector > > nullspace; - lr->getNullSpace(nullspace); + lrs->getNullSpace(nullspace); std::cout << "Nullspace: " << nullspace.size() << " x " << nullspace[0].size() << std::endl; cout << "Number of null vectors: " << nullspace.size() << endl; cout << "Vector sizes: " << nullspace[0].size() << endl; diff --git a/src/LRSplineVolume.cpp b/src/LRSplineVolume.cpp index 652e1f7..2c1d7f2 100644 --- a/src/LRSplineVolume.cpp +++ b/src/LRSplineVolume.cpp @@ -1591,6 +1591,29 @@ bool LRSplineVolume::isLinearIndepByOverloading(bool verbose) { iterationCount++; } while((uint) lastOverloadCount != overloaded.size()); + if(verbose && overloaded.size() > 0) { + std::cout << "The following ( " << overloaded.size() << ") basis functions could not be ruled out from linear dependence:" << std::endl; + double min[3] = { 1e99, 1e99, 1e99}; + double max[3] = {-1e99, -1e99, -1e99}; + std::set domain; + for(auto b : overloaded) { + std::cout << *b << std::endl; + for(auto el : b->support()) { + domain.insert(el->getId()); + } + for(int i=0; i<3; i++) { + min[i] = std::min(min[i], b->getParmin(i)); + max[i] = std::max(max[i], b->getParmax(i)); + } + } + std::cout << "These have support on the following (" << domain.size() << ") elements" << std::endl; + for(auto iel : domain) { + std::cout << *element_[iel] << std::endl; + } + std::cout << "These are contained in the bounding box:" << std::endl; + std::cout << " Lower Left: " << min[0] << ", " << min[1] << ", " << min[2] << std::endl; + std::cout << " Upper Right: " << max[0] << ", " << max[1] << ", " << max[2] << std::endl; + } return overloaded.size() == 0; } From c1612463e17a9b531d4fcfabc2365d33c0310e6f Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Sun, 17 Jan 2021 11:19:33 +0100 Subject: [PATCH 06/11] BUGFIX: Volume refinement was wrong (overlapping MeshRectangle algorithm) --- include/LRSpline/MeshRectangle.h | 2 +- src/LRSplineVolume.cpp | 154 ++++++++++++------- src/MeshRectangle.cpp | 244 +++++++++++++++---------------- 3 files changed, 221 insertions(+), 179 deletions(-) diff --git a/include/LRSpline/MeshRectangle.h b/include/LRSpline/MeshRectangle.h index 4bf454e..ae0863f 100644 --- a/include/LRSpline/MeshRectangle.h +++ b/include/LRSpline/MeshRectangle.h @@ -44,7 +44,6 @@ class MeshRectangle : public Streamable { bool contains(const MeshRectangle *rect) const; bool splits(Basisfunction *basis) const; bool splits(Element *el) const; - int makeOverlappingRects(std::vector &newGuys, int meshIndex, bool allowSplits) ; int multiplicity() const { return multiplicity_; }; int constDirection() const; @@ -56,6 +55,7 @@ class MeshRectangle : public Streamable { virtual void write(std::ostream &os) const; static bool addUniqueRect(std::vector &rects, MeshRectangle* newRect); + static int makeOverlappingRects(MeshRectangle* first, MeshRectangle* second, MeshRectangle **new1, MeshRectangle **new2); // private: std::vector start_; diff --git a/src/LRSplineVolume.cpp b/src/LRSplineVolume.cpp index 2c1d7f2..a575ab5 100644 --- a/src/LRSplineVolume.cpp +++ b/src/LRSplineVolume.cpp @@ -1227,6 +1227,7 @@ bool LRSplineVolume::enforceIsotropic() { } MeshRectangle* LRSplineVolume::insert_line(MeshRectangle *newRect) { + // Error test input if(newRect->start_[0] < start_[0] || newRect->start_[1] < start_[1] || newRect->start_[2] < start_[2] || @@ -1250,67 +1251,118 @@ MeshRectangle* LRSplineVolume::insert_line(MeshRectangle *newRect) { #ifdef TIME_LRSPLINE PROFILE("meshrectangle verification"); #endif - for(uint i=0; imakeOverlappingRects(newGuys, j, true); - if(status == 1) { // deleted j, i kept unchanged - j--; - continue; - } else if(status == 2) { // j kept unchanged, delete i - delete meshrect_[i]; - meshrect_.erase(meshrect_.begin() + i); - i--; - break; - } else if(status == 3) { // j kept unchanged, i added to newGuys - meshrect_.erase(meshrect_.begin() + i); - i--; - break; - } else if(status == 4) { // deleted j, i added to newGuys - meshrect_.erase(meshrect_.begin() + i); - i--; - j--; - break; - } else if(status == 5) { // deleted j, i duplicate in newGuys - delete meshrect_[i]; - meshrect_.erase(meshrect_.begin() + i); - i--; - j--; - break; - } else if(status == 6) { // j kept unchanged, i duplicate in newGuys - delete meshrect_[i]; - meshrect_.erase(meshrect_.begin() + i); - i--; - break; - } - } - } bool change = true; + bool did_del_j = false; + MeshRectangle *new1, *new2; while(change) { change = false; - for(uint i=0; ii; j++) { - int status = newGuys[i]->makeOverlappingRects(newGuys, j, false); - if(status == 1) { //deleted j, i kept unchanged - ; - } else if(status == 2) { // j kept unchanged, deleted i - delete newGuys[i]; - newGuys.erase(newGuys.begin() + i); - } else if(status == 3) { // j kept unchanged, i added to newGuys - newGuys.erase(newGuys.begin() + i); - } else if(status == 4) { // deleted j, i added to newGuys + for(uint j=0; !change && jequals(new1) || m->contains(new1)) + exists = true; + if(exists || !MeshRectangle::addUniqueRect(newGuys, new1)) + delete new1; + else + change = true; + } else if(status == 7) { // i,j unchanged, added both new1 and new2 + bool exists = false; + for(auto m : meshrect_) + if(m->equals(new1) || m->contains(new1)) + exists = true; + if(exists || !MeshRectangle::addUniqueRect(newGuys, new1)) + delete new1; + else + change = true; + exists = false; + for(auto m : meshrect_) + if(m->equals(new2) || m->contains(new2)) + exists = true; + if(exists || !MeshRectangle::addUniqueRect(newGuys, new2)) + delete new2; + else + change = true; + } + } + if(change) break; + if(did_del_j) continue; + + for(uint i=j+1; !change && i 0) { + i--; change = true; - break; + } else if(status == 4) { // j changed, i unchanged + change = true; + } else if(status == 5) { // j unchanged, i changed + change = true; + } else if(status == 6) { // j unchanged, i unchanged, new1 added + bool exists = false; + for(auto m : meshrect_) + if(m->equals(new1) || m->contains(new1)) + exists = true; + if(exists || !MeshRectangle::addUniqueRect(newGuys, new1)) + delete new1; + else + change = true; + } else if(status == 7) { // i,j unchanged, added both new1 and new2 + bool exists1 = false; + bool exists2 = false; + bool did_add_func = false; + for(auto m : meshrect_) { + if(m->equals(new1) || m->contains(new1)) + exists1 = true; + if(m->equals(new2) || m->contains(new2)) + exists2 = true; + } + if(exists1 || !MeshRectangle::addUniqueRect(newGuys, new1)) + delete new1; + else + change = true; + if(exists2 || !MeshRectangle::addUniqueRect(newGuys, new2)) + delete new2; + else + change = true; } } } } + } // end meshrectangle verification timer HashSet newFuncStp1, newFuncStp2; diff --git a/src/MeshRectangle.cpp b/src/MeshRectangle.cpp index 04dafa6..c0de60f 100644 --- a/src/MeshRectangle.cpp +++ b/src/MeshRectangle.cpp @@ -13,6 +13,8 @@ bool MeshRectangle::addUniqueRect(std::vector &rects, MeshRectan for(MeshRectangle *m : rects) { if(m->equals(newRect) ) { return false; + } else if(m->contains(newRect)) { + return false; } } rects.push_back(newRect); @@ -131,25 +133,22 @@ bool MeshRectangle::contains(const MeshRectangle *rect) const { return false; } -int MeshRectangle::makeOverlappingRects(std::vector &newGuys, int meshIndex, bool allowSplits) { - int c1 = constDir_; // constant index +int MeshRectangle::makeOverlappingRects(MeshRectangle* first, MeshRectangle* second, MeshRectangle **new1, MeshRectangle **new2) { + int c1 = first->constDir_; // constant index int v1 = (c1+1)%3; // first variable index int v2 = (c1+2)%3; // second variable index int v[] = {v1, v2}; // index way of referencing these bool addThisToNewGuys = false; - MeshRectangle *rect = newGuys.at(meshIndex); - if( ! this->overlaps(rect) ) + if( ! first->overlaps(second) ) return 0; - if( this->contains(rect) ) { - newGuys.erase(newGuys.begin() + meshIndex); - // std::cout << "Deleted: " << *rect << std::endl; - // std::cout << " contained in : " << *this << std::endl; - delete rect; + if( first->contains(second) ) { + // std::cout << "CONTAINED: " << *second << std::endl; + // std::cout << " contained in : " << *first << std::endl; return 1; } - if( rect->contains(this) ) { - // std::cout << "Deleted: " << *this << std::endl; - // std::cout << " contained in : " << *rect << std::endl; + if( second->contains(first) ) { + // std::cout << "CONTAINED: " << *first << std::endl; + // std::cout << " contained in : " << *second << std::endl; return 2; } @@ -167,94 +166,98 @@ int MeshRectangle::makeOverlappingRects(std::vector &newGuys, in * +-------+------+ */ // if the two mesh rectangles perfectly line up, keep only one of them - if((stop_[v[i]] == rect->start_[v[i]] || - start_[v[i]] == rect->stop_[v[i]] )) { - - if(start_[v[j]] == rect->start_[v[j]] && - stop_[v[j]] == rect->stop_[v[j]] ) { - double min = std::min(start_[v[i]], rect->start_[v[i]]); - double max = std::max(stop_[v[i]], rect->stop_[v[i]] ); - // std::cout << "Deleted: " << *rect << std::endl; - // std::cout << " merged with : " << *this << std::endl; - newGuys.erase(newGuys.begin() + meshIndex); - delete rect; - start_[v[i]] = min; - stop_[v[i]] = max; - if(addUniqueRect(newGuys, this)) - return 4; - else - return 5; - /* - * ADD ELONGATED RECT - * y3 +-------+ - * y2 | +------+ +-------+------+ - * | | | => | | - * y1 +-------+ | +-------+------+ - * | | new one - * y0 +------+ - * x0 x1 x3 - */ - - } else if((start_[v[j]] < rect->start_[v[j]] && - stop_[v[j]] < rect->stop_[v[j]] ) || - (start_[v[j]] > rect->start_[v[j]] && - stop_[v[j]] > rect->stop_[v[j]] )) { - double x0 = std::min(rect->start_[v[i]], start_[v[i]]); - double x3 = std::max(rect->stop_[v[i]], stop_[v[i]] ); - double y1 = std::max(rect->start_[v[j]], start_[v[j]]); - double y2 = std::min(rect->stop_[v[j]], stop_[v[j]] ); + if((first->stop_[v[i]] == second->start_[v[i]] || + first->start_[v[i]] == second->stop_[v[i]] )) { + + if(first->start_[v[j]] == second->start_[v[j]] && + first->stop_[v[j]] == second->stop_[v[j]] ) { + double min = std::min(first->start_[v[i]], second->start_[v[i]]); + double max = std::max(first->stop_[v[i]], second->stop_[v[i]] ); + // std::cout << "ELONGATED: " << *first << " (request delete)" << std::endl; + // std::cout << " merged (elongated) with : " << *second << " (kept)"; + second->start_[v[i]] = min; + second->stop_[v[i]] = max; + // std::cout << " new support: " << *second << std::endl; + return 3; + /* + * ADD ELONGATED RECT + * y3 +-------+ + * y2 | +------+ +-------+------+ + * | | | => | | + * y1 +-------+ | +-------+------+ + * | | new one + * y0 +------+ + * x0 x1 x3 + */ + + } else if((first->start_[v[j]] < second->start_[v[j]] && + first->stop_[v[j]] < second->stop_[v[j]] ) || + (first->start_[v[j]] > second->start_[v[j]] && + first->stop_[v[j]] > second->stop_[v[j]] )) { + double x0 = std::min(second->start_[v[i]], first->start_[v[i]]); + double x3 = std::max(second->stop_[v[i]], first->stop_[v[i]] ); + double y1 = std::max(second->start_[v[j]], first->start_[v[j]]); + double y2 = std::min(second->stop_[v[j]], first->stop_[v[j]] ); double start[3]; double stop[3]; - start[c1] = start_[c1]; - stop[c1] = stop_[c1]; + start[c1] = first->start_[c1]; + stop[c1] = first->stop_[c1]; start[v[i]] = x0; stop[v[i]] = x3; start[v[j]] = y1; stop[v[j]] = y2; - if(allowSplits) { - MeshRectangle *m1 = new MeshRectangle(start, stop); - if(!addUniqueRect(newGuys, m1)) - delete m1; - } + *new1 = new MeshRectangle(start, stop); + // std::cout << "Creating intersecting rect between " << *second << " and " << *first << " Result : " << **new1 << std::endl; + return 6; } } /* - * EXPAND 'RECT' (RIGHT ONE) + * EXPAND 'second' (RIGHT ONE) * +-------+ * | +--+------+ * | | | | * | +--+------+ * +-------+ */ - if(start_[v[i]] < rect->start_[v[i]] && - start_[v[j]] <= rect->start_[v[j]] && - stop_[v[j]] >= rect->stop_[v[j]]) { // expand the support of rect DOWN in v[i] - rect->start_[v[i]] = start_[v[i]]; + if(first->start_[v[i]] < second->start_[v[i]] && + first->start_[v[j]] <= second->start_[v[j]] && + first->stop_[v[j]] >= second->stop_[v[j]]) { // expand the support of second DOWN in v[i] + // std::cout << "Extended support (second-DOWN). Old Support: " << *second ; + second->start_[v[i]] = first->start_[v[i]]; + // std::cout << "new support: " << *second << std::endl; + return 4; } - if(stop_[v[i]] > rect->stop_[v[i]] && - start_[v[j]] <= rect->start_[v[j]] && - stop_[v[j]] >= rect->stop_[v[j]]) { // expand the support of rect UP in v[i] - rect->stop_[v[i]] = stop_[v[i]]; + if(first->stop_[v[i]] > second->stop_[v[i]] && + first->start_[v[j]] <= second->start_[v[j]] && + first->stop_[v[j]] >= second->stop_[v[j]]) { // expand the support of second UP in v[i] + // std::cout << "Extended support (second-UP). Old Support: " << *second ; + second->stop_[v[i]] = first->stop_[v[i]]; + // std::cout << "new support: " << *second << std::endl; + return 4; } /* - * EXPAND 'THIS' (LEFT ONE) + * EXPAND 'first' (LEFT ONE) * +------+ * +-----+--+ | * | | | | * +-----+--+ | * +------+ */ - if(rect->start_[v[i]] < start_[v[i]] && - rect->start_[v[j]] <= start_[v[j]] && - rect->stop_[v[j]] >= stop_[v[j]]) { - start_[v[i]] = rect->start_[v[i]]; - addThisToNewGuys = true; + if(second->start_[v[i]] < first->start_[v[i]] && + second->start_[v[j]] <= first->start_[v[j]] && + second->stop_[v[j]] >= first->stop_[v[j]]) { + // std::cout << "Extended support (first-DOWN). Old Support: " << *first ; + first->start_[v[i]] = second->start_[v[i]]; + // std::cout << "new support: " << *first << std::endl; + return 5; } - if(rect->stop_[v[i]] > stop_[v[i]] && - rect->start_[v[j]] <= start_[v[j]] && - rect->stop_[v[j]] >= stop_[v[j]]) { - stop_[v[i]] = rect->stop_[v[i]]; - addThisToNewGuys = true; + if(second->stop_[v[i]] > first->stop_[v[i]] && + second->start_[v[j]] <= first->start_[v[j]] && + second->stop_[v[j]] >= first->stop_[v[j]]) { + // std::cout << "Extended support (first-UP). Old Support: " << *first ; + first->stop_[v[i]] = second->stop_[v[i]]; + // std::cout << "new support: " << *first << std::endl; + return 5; } } @@ -270,14 +273,14 @@ int MeshRectangle::makeOverlappingRects(std::vector &newGuys, in * +--+ +--+ * note that this is after fixes above */ - if((rect->stop_[v1] >= stop_[v1] && - rect->start_[v1] <= start_[v1] && - stop_[v2] >= rect->stop_[v2] && - start_[v2] <= rect->start_[v2]) || - ( stop_[v1] >= rect->stop_[v1] && - start_[v1] <= rect->start_[v1] && - rect->stop_[v2] >= stop_[v2] && - rect->start_[v2] <= start_[v2])) { + if((second->stop_[v1] >= first->stop_[v1] && + second->start_[v1] <= first->start_[v1] && + first->stop_[v2] >= second->stop_[v2] && + first->start_[v2] <= second->start_[v2])|| + ( first->stop_[v1] >= second->stop_[v1] && + first->start_[v1] <= second->start_[v1] && + second->stop_[v2] >= first->stop_[v2] && + second->start_[v2]<= first->start_[v2])) { // ... ; /* @@ -292,53 +295,40 @@ int MeshRectangle::makeOverlappingRects(std::vector &newGuys, in * x0 x1 x2 x3 x1 x2 * these are the two new ones */ - } else { - double x0 = std::min(rect->start_[v1], start_[v1]); - double x1 = std::max(rect->start_[v1], start_[v1]); - double x2 = std::min(rect->stop_[v1], stop_[v1] ); - double x3 = std::max(rect->stop_[v1], stop_[v1] ); - double y0 = std::min(rect->start_[v2], start_[v2]); - double y1 = std::max(rect->start_[v2], start_[v2]); - double y2 = std::min(rect->stop_[v2], stop_[v2] ); - double y3 = std::max(rect->stop_[v2], stop_[v2] ); + } else if(first->start_[v1] < second->stop_[v1] && first->stop_[v1] > second->start_[v1] && + first->start_[v2] < second->stop_[v2] && first->stop_[v2] > second->start_[v2]) { + double x0 = std::min(second->start_[v1], first->start_[v1]); + double x1 = std::max(second->start_[v1], first->start_[v1]); + double x2 = std::min(second->stop_[v1], first->stop_[v1] ); + double x3 = std::max(second->stop_[v1], first->stop_[v1] ); + double y0 = std::min(second->start_[v2], first->start_[v2]); + double y1 = std::max(second->start_[v2], first->start_[v2]); + double y2 = std::min(second->stop_[v2], first->stop_[v2] ); + double y3 = std::max(second->stop_[v2], first->stop_[v2] ); double start[3]; double stop[3]; - start[c1] = start_[c1]; - stop[c1] = stop_[c1]; - - if(allowSplits) { - start[v1] = x0; - stop[v1] = x3; - start[v2] = y1; - stop[v2] = y2; - MeshRectangle *m1 = new MeshRectangle(start, stop); - - start[v1] = x1; - stop[v1] = x2; - start[v2] = y0; - stop[v2] = y3; - MeshRectangle *m2 = new MeshRectangle(start, stop); - - // std::cout << "Added: " << *m1 << std::endl; - // std::cout << "Added: " << *m2 << std::endl; - // std::cout << " overlaps from : " << *rect << std::endl; - // std::cout << " overlaps from : " << *this << std::endl; - - if(!addUniqueRect(newGuys, m1)) - delete m1; - if(!addUniqueRect(newGuys, m2)) - delete m2; - } - - } - if(addThisToNewGuys) { - // std::cout << "Moved: " << *this << std::endl; - if(addUniqueRect(newGuys, this)) - return 3; - else - return 6; + start[c1] = first->start_[c1]; + stop[c1] = first->stop_[c1]; + + start[v1] = x0; + stop[v1] = x3; + start[v2] = y1; + stop[v2] = y2; + *new1 = new MeshRectangle(start, stop); + + start[v1] = x1; + stop[v1] = x2; + start[v2] = y0; + stop[v2] = y3; + *new2 = new MeshRectangle(start, stop); + + // std::cout << "Added: " << **new1 << std::endl; + // std::cout << "Added: " << **new2 << std::endl; + // std::cout << " overlaps from : " << *second << std::endl; + // std::cout << " overlaps from : " << *first << std::endl; + return 7; } - return 0; + return -1; } bool MeshRectangle::splits(Element *el) const { From 01c1ff35a76b970ff2873cf938dbf47ac0dc43de Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Tue, 19 Jan 2021 15:51:06 +0100 Subject: [PATCH 07/11] Tests: added more tests for volume refinement and MeshRectangle behaviour --- test/RefinementUnchanged/extend_meshrect.inp | 3 +++ test/RefinementUnchanged/extend_meshrect.reg | 9 +++++++++ test/RefinementUnchanged/overlap_w.inp | 4 ++++ test/RefinementUnchanged/overlap_w.reg | 9 +++++++++ .../volume_all_coved_by_one_meshrect.inp | 12 ++++++++++++ .../volume_all_coved_by_one_meshrect.reg | 9 +++++++++ test/RefinementUnchanged/volume_merging_rect.inp | 4 ++++ test/RefinementUnchanged/volume_merging_rect.reg | 9 +++++++++ .../volume_two_rect_corner_overlap.inp | 3 +++ .../volume_two_rect_corner_overlap.reg | 9 +++++++++ 10 files changed, 71 insertions(+) create mode 100644 test/RefinementUnchanged/extend_meshrect.inp create mode 100644 test/RefinementUnchanged/extend_meshrect.reg create mode 100644 test/RefinementUnchanged/overlap_w.inp create mode 100644 test/RefinementUnchanged/overlap_w.reg create mode 100644 test/RefinementUnchanged/volume_all_coved_by_one_meshrect.inp create mode 100644 test/RefinementUnchanged/volume_all_coved_by_one_meshrect.reg create mode 100644 test/RefinementUnchanged/volume_merging_rect.inp create mode 100644 test/RefinementUnchanged/volume_merging_rect.reg create mode 100644 test/RefinementUnchanged/volume_two_rect_corner_overlap.inp create mode 100644 test/RefinementUnchanged/volume_two_rect_corner_overlap.reg diff --git a/test/RefinementUnchanged/extend_meshrect.inp b/test/RefinementUnchanged/extend_meshrect.inp new file mode 100644 index 0000000..4b9f392 --- /dev/null +++ b/test/RefinementUnchanged/extend_meshrect.inp @@ -0,0 +1,3 @@ +2 +0 1.5 0 1 1 4 1 +0 1.5 0 1 3 5 1 diff --git a/test/RefinementUnchanged/extend_meshrect.reg b/test/RefinementUnchanged/extend_meshrect.reg new file mode 100644 index 0000000..06e71c1 --- /dev/null +++ b/test/RefinementUnchanged/extend_meshrect.reg @@ -0,0 +1,9 @@ +-p1 3 -p2 3 -p3 3 -n1 10 -n2 4 -n3 8 extend_meshrect.inp + + all assertions passed + +Key LR-spline information: + number of mesh lines : 20 + number of elements : 100 + + diff --git a/test/RefinementUnchanged/overlap_w.inp b/test/RefinementUnchanged/overlap_w.inp new file mode 100644 index 0000000..6d4df63 --- /dev/null +++ b/test/RefinementUnchanged/overlap_w.inp @@ -0,0 +1,4 @@ +3 +2 1.5 1 4 0 3 1 +2 1.5 3 6 2 5 1 +2 1.5 0 7 1 4 1 diff --git a/test/RefinementUnchanged/overlap_w.reg b/test/RefinementUnchanged/overlap_w.reg new file mode 100644 index 0000000..83bf172 --- /dev/null +++ b/test/RefinementUnchanged/overlap_w.reg @@ -0,0 +1,9 @@ +-p1 3 -p2 3 -p3 3 -n1 10 -n2 8 -n3 4 overlap_w.inp + + all assertions passed + +Key LR-spline information: + number of basis functions: 331 + number of mesh lines : 23 + number of elements : 123 + diff --git a/test/RefinementUnchanged/volume_all_coved_by_one_meshrect.inp b/test/RefinementUnchanged/volume_all_coved_by_one_meshrect.inp new file mode 100644 index 0000000..65c24d7 --- /dev/null +++ b/test/RefinementUnchanged/volume_all_coved_by_one_meshrect.inp @@ -0,0 +1,12 @@ +11 +1 1.5 1 2 1 2 1 +1 1.5 3 4 1 2 1 +1 1.5 2 3 1 2 1 +1 1.5 0 4 0 3 1 +1 1.5 2 3 1 3 1 +1 1.5 3 4 2 3 1 +1 1.5 1 4 1 4 1 +1 1.5 3 5 2 5 1 +1 1.5 1 4 0 6 1 +1 1.5 0 5 0 6 1 +1 1.5 0 5 0 6 1 diff --git a/test/RefinementUnchanged/volume_all_coved_by_one_meshrect.reg b/test/RefinementUnchanged/volume_all_coved_by_one_meshrect.reg new file mode 100644 index 0000000..e708586 --- /dev/null +++ b/test/RefinementUnchanged/volume_all_coved_by_one_meshrect.reg @@ -0,0 +1,9 @@ +-p1 3 -p2 3 -p3 3 -n1 10 -n2 4 -n3 8 volume_all_coved_by_one_meshrect.inp + + + all assertions passed + +Key LR-spline information: + number of mesh lines : 20 + number of elements : 126 + diff --git a/test/RefinementUnchanged/volume_merging_rect.inp b/test/RefinementUnchanged/volume_merging_rect.inp new file mode 100644 index 0000000..d91439d --- /dev/null +++ b/test/RefinementUnchanged/volume_merging_rect.inp @@ -0,0 +1,4 @@ +3 +1 1.5 1 2 1 2 1 +1 1.5 3 4 1 2 1 +1 1.5 2 3 1 2 1 diff --git a/test/RefinementUnchanged/volume_merging_rect.reg b/test/RefinementUnchanged/volume_merging_rect.reg new file mode 100644 index 0000000..ad75e4b --- /dev/null +++ b/test/RefinementUnchanged/volume_merging_rect.reg @@ -0,0 +1,9 @@ +-p1 3 -p2 3 -p3 3 -n1 10 -n2 4 -n3 8 volume_merging_rect.inp + + all assertions passed + +Key LR-spline information: + number of basis functions: 320 + number of mesh lines : 20 + number of elements : 99 + diff --git a/test/RefinementUnchanged/volume_two_rect_corner_overlap.inp b/test/RefinementUnchanged/volume_two_rect_corner_overlap.inp new file mode 100644 index 0000000..e4ddb91 --- /dev/null +++ b/test/RefinementUnchanged/volume_two_rect_corner_overlap.inp @@ -0,0 +1,3 @@ +4 +1 1.5 1 4 2 5 1 +1 1.5 2 5 1 4 1 diff --git a/test/RefinementUnchanged/volume_two_rect_corner_overlap.reg b/test/RefinementUnchanged/volume_two_rect_corner_overlap.reg new file mode 100644 index 0000000..ef1167a --- /dev/null +++ b/test/RefinementUnchanged/volume_two_rect_corner_overlap.reg @@ -0,0 +1,9 @@ +-p1 3 -p2 3 -p3 3 -n1 10 -n2 4 -n3 8 volume_two_rect_corner_overlap.inp + + all assertions passed + +Key LR-spline information: + number of basis functions: 322 + number of mesh lines : 23 + number of elements : 110 + From cba6956dccff3e57d9a90f7f1d518e54fabe5fbb Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Tue, 19 Jan 2021 15:52:04 +0100 Subject: [PATCH 08/11] bump version 1.7.0 --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a846dcb..8c20100 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,8 +2,8 @@ PROJECT(LRSpline) CMAKE_MINIMUM_REQUIRED(VERSION 2.6) SET(LRSpline_VERSION_MAJOR 1) -SET(LRSpline_VERSION_MINOR 6) -SET(LRSpline_VERSION_PATCH 1) +SET(LRSpline_VERSION_MINOR 7) +SET(LRSpline_VERSION_PATCH 0) SET(LRSpline_VERSION ${LRSpline_VERSION_MAJOR}.${LRSpline_VERSION_MINOR}.${LRSpline_VERSION_PATCH}) IF(NOT WIN32) SET(TIME_LRSPLINE 1) From ba55872cb47ac8c6bc110d9ec077f6ce856b3665 Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Thu, 11 Jun 2020 15:44:42 +0200 Subject: [PATCH 09/11] Fixes #39: adds control over Element enumeration for consistency tests --- Apps/TestReadWrite.cpp | 3 +++ include/LRSpline/Element.h | 2 ++ include/LRSpline/LRSpline.h | 2 ++ src/Element.cpp | 12 ++++++++++++ src/LRSpline.cpp | 8 ++++++++ 5 files changed, 27 insertions(+) diff --git a/Apps/TestReadWrite.cpp b/Apps/TestReadWrite.cpp index 469eb54..35a3e5f 100644 --- a/Apps/TestReadWrite.cpp +++ b/Apps/TestReadWrite.cpp @@ -132,6 +132,9 @@ int main(int argc, char **argv) { inputfile >> *lrs; } + if(vol) lrv->renumberElements(); + else lrs->renumberElements(); + // test writing to file ofstream lrfile; lrfile.open("TestReadWrite.lr"); diff --git a/include/LRSpline/Element.h b/include/LRSpline/Element.h index 8f070da..dec35fd 100644 --- a/include/LRSpline/Element.h +++ b/include/LRSpline/Element.h @@ -98,6 +98,8 @@ class Element : public Streamable { virtual void read(std::istream &is); virtual void write(std::ostream &os) const; + bool operator<(const Element &other) const; + private: std::vector min; // lower left corner in typical 2 or 3 dimensions std::vector max; // upper right corner diff --git a/include/LRSpline/LRSpline.h b/include/LRSpline/LRSpline.h index f1d258a..099e645 100644 --- a/include/LRSpline/LRSpline.h +++ b/include/LRSpline/LRSpline.h @@ -47,6 +47,8 @@ class LRSpline : public Streamable { virtual void generateIDs() const; + virtual void renumberElements() ; + // common get methods //! \brief returns the number of B-splines (basisfunctions) in this LR-spline object diff --git a/src/Element.cpp b/src/Element.cpp index d4dde5e..0a3c4b8 100644 --- a/src/Element.cpp +++ b/src/Element.cpp @@ -298,6 +298,18 @@ bool Element::isOverloaded() const { } } return false; + +} + +bool Element::operator<(const Element &other) const { + for(size_t i=0; imin[i] > other.getParmin(i)) + return false; + else if(this->min[i] < other.getParmin(i)) + return true; + // else(this->min[i] == other.getParmin(i)) continue + } + return false; } } // end namespace LR diff --git a/src/LRSpline.cpp b/src/LRSpline.cpp index 9794490..20840d8 100644 --- a/src/LRSpline.cpp +++ b/src/LRSpline.cpp @@ -11,6 +11,14 @@ LRSpline::LRSpline() { element_.resize(0); } +bool compare(const Element* a, const Element* b) {return *a < *b;}; + +void LRSpline::renumberElements() { + std::sort(element_.begin(), element_.end(), compare); + for(size_t i=0; isetId(i); +} + void LRSpline::generateIDs() const { uint i=0; for(Basisfunction *b : basis_) From 5bb939f1ecb547b6f5d1e3d57ea90261c9d04049 Mon Sep 17 00:00:00 2001 From: Kjetil Johannessen Date: Thu, 11 Jun 2020 16:28:28 +0200 Subject: [PATCH 10/11] Fixes #40: orders basisfunction index by sorting on the smallest differing knot from the local knot vectors --- Apps/TestReadWrite.cpp | 3 +++ include/LRSpline/Basisfunction.h | 1 + include/LRSpline/LRSpline.h | 3 ++- src/Basisfunction.cpp | 15 ++++++++++++ src/LRSpline.cpp | 41 +++++++++++++++++++++++++++++--- 5 files changed, 59 insertions(+), 4 deletions(-) diff --git a/Apps/TestReadWrite.cpp b/Apps/TestReadWrite.cpp index 35a3e5f..1507c04 100644 --- a/Apps/TestReadWrite.cpp +++ b/Apps/TestReadWrite.cpp @@ -135,6 +135,9 @@ int main(int argc, char **argv) { if(vol) lrv->renumberElements(); else lrs->renumberElements(); + if(vol) lrv->renumberBasisfunctions(); + else lrs->renumberBasisfunctions(); + // test writing to file ofstream lrfile; lrfile.open("TestReadWrite.lr"); diff --git a/include/LRSpline/Basisfunction.h b/include/LRSpline/Basisfunction.h index 2e42097..c0ce1f5 100644 --- a/include/LRSpline/Basisfunction.h +++ b/include/LRSpline/Basisfunction.h @@ -179,6 +179,7 @@ class Basisfunction : public Streamable { // operator overloading bool equals(const Basisfunction &other) const ; bool operator==(const Basisfunction &other) const; + bool operator<( const Basisfunction &other) const; void operator+=(const Basisfunction &other) ; std::vector& operator[](int i) { return knots_[i]; } ; const std::vector& operator[](int i) const { return knots_[i]; } ; diff --git a/include/LRSpline/LRSpline.h b/include/LRSpline/LRSpline.h index 099e645..4914fcc 100644 --- a/include/LRSpline/LRSpline.h +++ b/include/LRSpline/LRSpline.h @@ -48,6 +48,7 @@ class LRSpline : public Streamable { virtual void generateIDs() const; virtual void renumberElements() ; + virtual void renumberBasisfunctions() ; // common get methods @@ -93,7 +94,7 @@ class LRSpline : public Streamable { // traditional get methods Element* getElement(int i) { return element_[i]; }; const Element* getElement(int i) const { return element_[i]; }; - Basisfunction* getBasisfunction(int iBasis) ; + Basisfunction* getBasisfunction(int iBasis); const Basisfunction* getBasisfunction(int iBasis) const ; // refinement functions diff --git a/src/Basisfunction.cpp b/src/Basisfunction.cpp index 66f89d5..8298ac3 100644 --- a/src/Basisfunction.cpp +++ b/src/Basisfunction.cpp @@ -695,6 +695,21 @@ bool Basisfunction::operator==(const Basisfunction &other) const { return equals(other); } +/************************************************************************************************************************//** + * \brief Test for B-spline ordering according to geometric factor instead of hash function + * \param other The other B-spline to check against + * \returns True if this.knot[i][j] < other.knot[i][j] for the smallest possible (i,j) whith different values + ***************************************************************************************************************************/ +bool Basisfunction::operator<(const Basisfunction &other) const { + for(uint i=0; i other[i][j]) return false; + else if(knots_[i][j] < other[i][j]) return true; + // else continue; + + return false; // all knot vectors identical +} + /************************************************************************************************************************//** * \brief Test for B-spline equality * \param other The other B-spline to check against diff --git a/src/LRSpline.cpp b/src/LRSpline.cpp index 20840d8..ba76f1b 100644 --- a/src/LRSpline.cpp +++ b/src/LRSpline.cpp @@ -11,19 +11,40 @@ LRSpline::LRSpline() { element_.resize(0); } -bool compare(const Element* a, const Element* b) {return *a < *b;}; +bool compareEl(const Element* a, const Element* b) {return *a < *b;}; void LRSpline::renumberElements() { - std::sort(element_.begin(), element_.end(), compare); + this->generateIDs(); + std::sort(element_.begin(), element_.end(), compareEl); for(size_t i=0; isetId(i); } +bool compareFunc(const Basisfunction* a, const Basisfunction* b) {return *a < *b;}; + +void LRSpline::renumberBasisfunctions() { + std::vector tmp(basis_.size()); + std::copy(basis_.begin(), basis_.end(), tmp.begin()); + std::sort(tmp.begin(), tmp.end(), compareFunc); + for(size_t i=0; isetId(i); +} + void LRSpline::generateIDs() const { + /* enumeration by hash function (old way), sensitive to scaling and bit precision uint i=0; for(Basisfunction *b : basis_) b->setId(i++); - for(i=0; i tmp(basis_.size()); + std::copy(basis_.begin(), basis_.end(), tmp.begin()); + std::sort(tmp.begin(), tmp.end(), compareFunc); + for(size_t i=0; isetId(i); + + for(size_t i=0; isetId(i); } @@ -123,5 +144,19 @@ void LRSpline::rebuildDimension(int dimvalue) { dim_ = dimvalue; } +Basisfunction* LRSpline::getBasisfunction(int iBasis) { + if(iBasis<0 || iBasis>=basis_.size()) return NULL; + for(auto b : basis_) + if(b->getId() == iBasis) return b; + return NULL; +} + +const Basisfunction* LRSpline::getBasisfunction(int iBasis) const { + if(iBasis<0 || iBasis>=basis_.size()) return NULL; + for(auto b : basis_) + if(b->getId() == iBasis) return b; + return NULL; +} + } // end namespace LR From e7a8635843ca39bce5156c8ce52daca9ad854c45 Mon Sep 17 00:00:00 2001 From: Eivind Fonn Date: Wed, 23 Sep 2020 17:49:24 +0200 Subject: [PATCH 11/11] Fix read/copy basis function numbering After basis function numbering changed, files no longer necessarily have the basis functions listed in order of ID, therefore we can't just update pointers based on file order and expect to get the right ones. --- src/LRSplineSurface.cpp | 6 +++--- src/LRSplineVolume.cpp | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/LRSplineSurface.cpp b/src/LRSplineSurface.cpp index 5cf914a..d23493b 100644 --- a/src/LRSplineSurface.cpp +++ b/src/LRSplineSurface.cpp @@ -207,7 +207,7 @@ LRSplineSurface::~LRSplineSurface() { LRSplineSurface* LRSplineSurface::copy() const { generateIDs(); - std::vector basisVector; + std::vector basisVector(nBasisFunctions()); // flat list to make it quicker to update pointers from Basisfunction to Element and back again LRSplineSurface *returnvalue = new LR::LRSplineSurface(); @@ -215,7 +215,7 @@ LRSplineSurface* LRSplineSurface::copy() const { for(Basisfunction* b : basis_) { Basisfunction *newB = b->copy(); returnvalue -> basis_.insert(newB); - basisVector.push_back(newB); + basisVector[b->getId()] = newB; } for(Element *e : element_) { @@ -2869,7 +2869,7 @@ void LRSplineSurface::read(std::istream &is) { Basisfunction *b = new Basisfunction(dim_, order_[0], order_[1]); b->read(is); basis_.insert(b); - basisVector[i] = b; + basisVector[b->getId()] = b; } // get rid of more comments and spaces diff --git a/src/LRSplineVolume.cpp b/src/LRSplineVolume.cpp index a575ab5..0df7f46 100644 --- a/src/LRSplineVolume.cpp +++ b/src/LRSplineVolume.cpp @@ -182,7 +182,7 @@ void LRSplineVolume::initMeta() { LRSplineVolume* LRSplineVolume::copy() const { generateIDs(); - std::vector basisVector; + std::vector basisVector(nBasisFunctions()); // flat list to make it quicker to update pointers from Basisfunction to Element and back again LRSplineVolume *returnvalue = new LR::LRSplineVolume(); @@ -190,7 +190,7 @@ LRSplineVolume* LRSplineVolume::copy() const { for(Basisfunction* b : basis_) { Basisfunction *newB = b->copy(); returnvalue -> basis_.insert(newB); - basisVector.push_back(newB); + basisVector[b->getId()] = newB; } for(Element *e : element_) { @@ -2242,7 +2242,7 @@ void LRSplineVolume::read(std::istream &is) { Basisfunction *b = new Basisfunction(dim_, 3, allOrder); b->read(is); basis_.insert(b); - basisVector[i] = b; + basisVector[b->getId()] = b; } // get rid of more comments and spaces