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/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..8c20100 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ PROJECT(LRSpline) CMAKE_MINIMUM_REQUIRED(VERSION 2.6) SET(LRSpline_VERSION_MAJOR 1) -SET(LRSpline_VERSION_MINOR 6) +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) @@ -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/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/include/LRSpline/Element.h b/include/LRSpline/Element.h index 69f547f..8f070da 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,18 @@ 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]; } + + // topological operators + bool touches(const Element* el) const ; + std::vector neighbours() ; void updateBasisPointers(std::vector &basis) ; @@ -96,7 +106,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/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/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/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/Element.cpp b/src/Element.cpp index c23a516..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; @@ -12,8 +13,9 @@ namespace LR { * \brief Default constructor ***************************************************************************************************************************/ Element::Element() { - id_ = -1; - overloadCount = 0; + id_ = -1; + overloadCount_ = 0; + level_ = std::vector(0); } /************************************************************************************************************************//** @@ -22,7 +24,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 +46,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 +75,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 +102,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 @@ -117,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 ***************************************************************************************************************************/ @@ -156,6 +202,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 +234,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 +270,9 @@ void Element::write(std::ostream &os) const { os << b->getId(); isFirst = false; } - os << "}"; + os << "} lvl " << level_[0]; + for(uint i=1; icloseGaps(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 @@ -1187,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] || @@ -1210,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; @@ -1551,6 +1643,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; } diff --git a/src/MeshRectangle.cpp b/src/MeshRectangle.cpp index 52223f8..c0de60f 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 { @@ -12,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); @@ -130,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; } @@ -166,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; } } @@ -269,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])) { // ... ; /* @@ -291,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 { 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 + diff --git a/test/TestLevelCounting/testSurf.reg b/test/TestLevelCounting/testSurf.reg new file mode 100644 index 0000000..d0cd828 --- /dev/null +++ b/test/TestLevelCounting/testSurf.reg @@ -0,0 +1,9 @@ +-p1 3 -p2 3 -n1 9 -n2 9 + +Length=0.0625 level 4, 4 (expects lvl 4) - OK +Length=0.125 level 3, 3 (expects lvl 3) - OK +Length=0.25 level 2, 2 (expects lvl 2) - OK +Length=0.5 level 1, 1 (expects lvl 1) - OK +Length=1 level 0, 0 (expects lvl 0) - OK +==================================================== +All assertions passed diff --git a/test/TestLevelCounting/testVol.reg b/test/TestLevelCounting/testVol.reg new file mode 100644 index 0000000..9127d8d --- /dev/null +++ b/test/TestLevelCounting/testVol.reg @@ -0,0 +1,8 @@ +-p1 3 -p2 3 -p3 3 -n1 9 -n2 9 -n3 9 -vol +Length=0.0625 level 4, 4, 4 (expects lvl 4) - OK +Length=0.125 level 3, 3, 3 (expects lvl 3) - OK +Length=0.25 level 2, 2, 2 (expects lvl 2) - OK +Length=0.5 level 1, 1, 1 (expects lvl 1) - OK +Length=1 level 0, 0, 0 (expects lvl 0) - OK +==================================================== +All assertions passed