Added force max T-joints for volumetric mesh

VikingScientist · Dec 22, 2020 · 8afbd4f · 8afbd4f
1 parent bb75482
commit 8afbd4f
Show file tree

Hide file tree

Showing 4 changed files with 85 additions and 2 deletions.
diff --git 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<Element*> neighbours() ;
+
 	void updateBasisPointers(std::vector<Basisfunction*> &basis) ;
 
 	virtual void read(std::istream &is);

diff --git 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<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

diff --git a/src/Element.cpp b/src/Element.cpp
@@ -3,6 +3,7 @@
 #include "LRSpline/Meshline.h"
 #include "LRSpline/Basisfunction.h"
 #include <stdlib.h>
+#include <set>
 
 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; dim<min.size(); dim++) {
+		if(el->getParmin(dim) == this->getParmax(dim) || el->getParmax(dim) == this->getParmin(dim)) {
+			bool overlaps = true;
+			for(uint i=0; i<dim; i++) {
+				if(i==dim) continue;
+				overlaps &= el->getParmax(i) > this->getParmin(i);
+				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*> Element::neighbours() {
+	std::vector<Element*> 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<Element*> 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
@@ -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<MeshRectangle*> 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