Max T joints

Apps/LinearIndep.cpp

@@ -5,6 +5,7 @@
 #include <fstream>
 #include <cstdlib>
 #include "LRSpline/LRSplineSurface.h"
+#include "LRSpline/LRSplineVolume.h"
 #include "LRSpline/Profiler.h"
 #include "LRSpline/Element.h"
 #include "LRSpline/Meshline.h"
@@ -92,19 +93,30 @@ int main(int argc, char **argv) {
 		exit(3);
 	}
 
-	LRSplineSurface *lr;
+	LRSpline        *lr;
+	LRSplineSurface *lrs;
+	LRSplineVolume  *lrv;
 
 	ifstream inputfile;
 	inputfile.open(inputFileName);
 	if(!inputfile.is_open()) {
 		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<Meshline*> *newLines = new vector<Meshline*>();
-		// 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; i<newLines->size(); 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<vector<boost::rational<long long> > > 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;

Apps/TestLevelCounting.cpp

@@ -0,0 +1,141 @@
+#include <stdio.h>
+#include <cstdlib>
+#include <iostream>
+#include <string.h>
+#include <fstream>
+#include <cmath>
+#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  <n>  polynomial ORDER (degree+1) in first parametric direction\n" \
+	                  "   -p2  <n>  polynomial order in second parametric direction\n" \
+	                  "   -p3  <n>  polynomial order in third parametric direction\n" \
+	                  "   -n1  <n>  number of basis functions in first parametric direction\n" \
+	                  "   -n2  <n>  number of basis functions in second parametric direction\n" \
+	                  "   -n3  <n>  number of basis functions in third parametric direction\n" \
+	                  "   -dim <n>  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<argc; i++) {
+		if(strcmp(argv[i], "-p1") == 0)
+			p1 = atoi(argv[++i]);
+		else if(strcmp(argv[i], "-p2") == 0)
+			p2 = atoi(argv[++i]);
+		else if(strcmp(argv[i], "-p3") == 0) {
+			p3 = atoi(argv[++i]);
+			vol = true;
+		} else if(strcmp(argv[i], "-n1") == 0)
+			n1 = atoi(argv[++i]);
+		else if(strcmp(argv[i], "-n2") == 0)
+			n2 = atoi(argv[++i]);
+		else if(strcmp(argv[i], "-n3") == 0) {
+			n3 = atoi(argv[++i]);
+			vol = true;
+		} else if(strcmp(argv[i], "-dim") == 0)
+			dim = atoi(argv[++i]);
+		else if(strcmp(argv[i], "-vol") == 0)
+			vol = true;
+		else if(strcmp(argv[i], "-help") == 0) {
+			cerr << "usage: " << argv[0] << " [parameters]" << endl << parameters.c_str();
+			exit(0);
+		} else {
+			cerr << "usage: " << argv[0] << " [parameters]" << endl << parameters.c_str();
+			exit(1);
+		}
+	}
+
+	// do some error testing on input
+	if(n1 < p1) {
+		cerr << "ERROR: n1 must be greater or equal to p1\n";
+		exit(2);
+	} else if(n2 < p2) {
+		cerr << "ERROR: n2 must be greater or equal to p2\n";
+		exit(2);
+	} else if(n3 < p3) {
+		cerr << "ERROR: n3 must be greater or equal to p3\n";
+		exit(2);
+	}
+
+
+	// read the surface, or make one up if none specified
+	LRSpline        *lr  = NULL;
+	LRSplineSurface *lrs = NULL;
+	LRSplineVolume  *lrv = NULL;
+	// make a uniform integer knot vector
+	std::vector<double> knot_u(n1 + p1);
+	std::vector<double> knot_v(n2 + p2);
+	std::vector<double> knot_w(n3 + p3);
+	for(int i=0; i<p1+n1; i++)
+		knot_u[i] = (i<p1) ? 0 : (i>n1) ? n1-p1+1 : i-p1+1;
+	for(int i=0; i<p2+n2; i++)
+		knot_v[i] = (i<p2) ? 0 : (i>n2) ? n2-p2+1 : i-p2+1;
+	for(int i=0; i<p3+n3; i++)
+		knot_w[i] = (i<p3) ? 0 : (i>n3) ? 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<double> cp(nCP);
+	int k=0;
+	for(int i=0; i<nCP; i++) // 839 as a generator over Z_853 gives a period of 425. Should suffice
+		cp[k++] = (i*839 % 853) / 853.0 + 0.1;  // rational weights also random and thus we need >0
+
+	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<double> corner(lr->nVariate(), 0.00001);
+		Element* el = lr->getElement(lr->getElementContaining(corner));
+		vector<int> 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; i<lr->nVariate(); i++) cout << ", " << el->getLevel(i);
+		cout << " (expects lvl " << expectedLevel << ")";
+		bool ok = true;
+		for(int i=0; i<lr->nVariate(); 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;
+}
+

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}")

include/LRSpline/Basisfunction.h

@@ -170,7 +170,7 @@ class Basisfunction : public Streamable {
 	std::vector<double>::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 ;
 

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<int>(dim, 0);
 	}
 	Element(std::vector<double> &lowerLeft, std::vector<double> &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<Element*> neighbours() ;
 
 	void updateBasisPointers(std::vector<Basisfunction*> &basis) ;
 
@@ -96,7 +106,10 @@ class Element : public Streamable {
 	HashSet<Basisfunction*> support_;
 	std::vector<int> 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<int> level_ ;
 
 };
 

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<MeshRectangle*>* newLines=NULL);
-	void enforceMaxTjoints(    std::vector<MeshRectangle*>* newLines=NULL);
 	void enforceMaxAspectRatio(std::vector<MeshRectangle*>* newLines=NULL);
 	*/
+	void aPosterioriFixes() ;
+	void enforceMaxTjoints();
 	bool enforceIsotropic();
 
 	// linear independence methods

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<MeshRectangle*> &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<MeshRectangle*> &rects, MeshRectangle* newRect);
+	static int makeOverlappingRects(MeshRectangle* first, MeshRectangle* second, MeshRectangle **new1, MeshRectangle **new2);
 
 // private:
 	std::vector<double> start_;