Speedup: LR evaluation by bezier extraction (also 3D boundary conditions)

src/ASM/ASMunstruct.h

@@ -115,7 +115,7 @@ public:
   virtual LR::LRSpline* evalSolution(const IntegrandBase& integrand) const = 0;
 
   //! \brief Returns a Bezier basis of order \a p.
-  static Go::BsplineBasis getBezierBasis(int p);
+  static Go::BsplineBasis getBezierBasis(int p, double start=-1, double end=1);
 
   //! \brief Returns a list of basis functions having support on given elements.
   IntVec getFunctionsForElements(const IntVec& elements);

src/ASM/LR/ASMLRSpline.C

@@ -230,11 +230,11 @@ bool ASMunstruct::refine (const LR::RefineData& prm,
 }
 
 
-Go::BsplineBasis ASMunstruct::getBezierBasis (int p)
+Go::BsplineBasis ASMunstruct::getBezierBasis (int p, double start, double end)
 {
   double knot[2*p];
-  std::fill(knot,   knot+p,  -1.0);
+  std::fill(knot,   knot+p,   start);
-  std::fill(knot+p, knot+2*p, 1.0);
+  std::fill(knot+p, knot+2*p, end);
   return Go::BsplineBasis(p,p,knot);
 }
 

src/ASM/LR/ASMu2D.C

@@ -861,21 +861,20 @@ size_t ASMu2D::getNoBoundaryElms (char lIndex, char ldim) const
   else if (ldim < 1 && lIndex > 0)
     return 1;
 
-  // Maybe there is a more elegent way to count the boundary elements, Kjetil?
+  LR::parameterEdge edge;
-  size_t nbel = 0;
+  switch(lIndex)
-  std::vector<LR::Element*>::iterator el = lrspline->elementBegin();
-  for (; el != lrspline->elementEnd(); el++)
   {
-    switch (lIndex)
+  case 1: edge = LR::WEST;  break;
-    {
+  case 2: edge = LR::EAST;  break;
-    case 1: if ((*el)->umin() == lrspline->startparam(0)) nbel++; break;
+  case 3: edge = LR::SOUTH; break;
-    case 2: if ((*el)->umax() == lrspline->endparam(0))   nbel++; break;
+  case 4: edge = LR::NORTH; break;
-    case 3: if ((*el)->vmin() == lrspline->startparam(1)) nbel++; break;
+  default:edge = LR::NONE;
-    case 4: if ((*el)->vmax() == lrspline->endparam(1))   nbel++; break;
-    }
   }
 
-  return nbel;
+  std::vector<LR::Element*> edgeElms;
+  lrspline->getEdgeElements(edgeElms, edge);
+  
+	return edgeElms.size();
 }
 
 

src/ASM/LR/ASMu3D.C

@@ -618,19 +618,19 @@ void ASMu3D::constrainLine (int fdir, int ldir, double xi, int dof,
 
 void ASMu3D::constrainCorner (int I, int J, int K, int dof, int code, char)
 {
-	std::cerr << "ASMu3D::constrainCorner not implemented properly yet" << std::endl;
+  int corner = LR::NONE;
-	exit(776654);
+  if(  I < 0) corner |= LR::WEST;
-#if 0
+  else        corner |= LR::EAST;
-	int n1, n2, n3;
+  if(  J < 0) corner |= LR::SOUTH;
-	if (!this->getSize(n1,n2,n3,1)) return;
+  else        corner |= LR::NORTH;
+  if(  K < 0) corner |= LR::BOTTOM;
+  else        corner |= LR::TOP;
 
-	int node = 1;
+  std::vector<LR::Basisfunction*> one_function;
-	if (I > 0) node += n1-1;
+  lrspline->getEdgeFunctions(one_function, LR::parameterEdge(corner));
-	if (J > 0) node += n1*(n2-1);
+  int node = one_function[0]->getId() + 1;
-	if (K > 0) node += n1*n2*(n3-1);
 
-	this->prescribe(node,dof,code);
+  this->prescribe(node, dof, code);
-#endif
 }
 
 
@@ -744,36 +744,23 @@ bool ASMu3D::updateCoords (const Vector& displ)
 size_t ASMu3D::getNoBoundaryElms (char lIndex, char ldim) const
 {
 	if (!lrspline) return 0;
-	/*TODO fix later (see ASMu2D::getNoBoundaryElms)
-	if (ldim < 1 && lIndex > 0)
-		return 1;
 
-	std::cerr <<"ASMu3D::getNoBoundaryElms not implemented properly yet" << std::endl;
+	LR::parameterEdge edge;
-	exit(776654);
+	switch(lIndex)
+	{
+	case 1: edge = LR::WEST;   break;
+	case 2: edge = LR::EAST;   break;
+	case 3: edge = LR::SOUTH;  break;
+	case 4: edge = LR::NORTH;  break;
+	case 5: edge = LR::BOTTOM; break;
+	case 6: edge = LR::TOP;    break;
+	default:edge = LR::NONE;
+	}
 
-#if  0
+	std::vector<LR::Element*> edgeElms;
-	else if (ldim < 2 && lIndex > 0 && lIndex <= 12)
+	lrspline->getEdgeElements(edgeElms, edge);
-		return lrspline->numCoefs((lIndex-1)/4) - lrspline->order((lIndex-1)/4) + 1;
 
-	int n1 = lrspline->numCoefs(0) - lrspline->order(0) + 1;
+	return edgeElms.size();
-	int n2 = lrspline->numCoefs(1) - lrspline->order(1) + 1;
-	int n3 = lrspline->numCoefs(2) - lrspline->order(2) + 1;
-
-	switch (lIndex)
-		{
-		case 1:
-		case 2:
-			return n2*n3;
-		case 3:
-		case 4:
-			return n1*n3;
-		case 5:
-		case 6:
-			return n1*n2;
-		}
-#endif
-	*/
-	return 0;
 }
 
 
@@ -915,9 +902,15 @@ bool ASMu3D::integrate (Integrand& integrand,
 	int p1 = lrspline->order(0);
 	int p2 = lrspline->order(1);
 	int p3 = lrspline->order(2);
-	Go::BsplineBasis basis1 = getBezierBasis(p1);
+	double u0 = lrspline->startparam(0);
-	Go::BsplineBasis basis2 = getBezierBasis(p2);
+	double u1 = lrspline->endparam(0);
-	Go::BsplineBasis basis3 = getBezierBasis(p3);
+	double v0 = lrspline->startparam(1);
+	double v1 = lrspline->endparam(1);
+	double w0 = lrspline->startparam(2);
+	double w1 = lrspline->endparam(2);
+	Go::BsplineBasis basis1 = getBezierBasis(p1, u0, u1);
+	Go::BsplineBasis basis2 = getBezierBasis(p2, v0, v1);
+	Go::BsplineBasis basis3 = getBezierBasis(p3, w0, w1);
 
 	Matrix BN(   p1*p2*p3, nGauss*nGauss*nGauss);
 	Matrix BdNdu(p1*p2*p3, nGauss*nGauss*nGauss);
@@ -930,9 +923,12 @@ bool ASMu3D::integrate (Integrand& integrand,
 				double u[2*p1];
 				double v[2*p2];
 				double w[2*p3];
-				basis1.computeBasisValues(xg[xi],   u, 1); 
+				double evalU = (xg[xi]  +1)*0.5*(u1-u0) + u0;
-				basis2.computeBasisValues(xg[eta],  v, 1);
+				double evalV = (xg[eta] +1)*0.5*(v1-v0) + v0;
-				basis3.computeBasisValues(xg[zeta], w, 1);
+				double evalW = (xg[zeta]+1)*0.5*(w1-w0) + w0;
+				basis1.computeBasisValues(evalU, u, 1); 
+				basis2.computeBasisValues(evalV, v, 1);
+				basis3.computeBasisValues(evalW, w, 1);
 				int ib=1; // basis function iterator
 				double sum = 0;
 				for(int k=0; k<p3; k++) {
@@ -1131,6 +1127,7 @@ bool ASMu3D::integrate (Integrand& integrand,
 					if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
 						evaluateBasis(fe, dNdu, d2Ndu2);
 					else
+						// evaluateBasis(fe, dNdu);
 						evaluateBasis(fe, dNdu, C, B) ;
 
 					// look for errors in bezier extraction
@@ -1234,6 +1231,8 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
 	const int t1 = 1 + abs(faceDir)%3; // first tangent direction
 	const int t2 = 1 + t1%3;           // second tangent direction
 
+  std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex);
+  size_t firstp = iit == firstBp.end() ? 0 : iit->second;
 	
 	int edge = LR::NONE;
 	if(lIndex == 1) 
@@ -1318,9 +1317,10 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
 
 		// --- Integration loop over all Gauss points in each direction --------
 
+		fe.iGP = firstp; // Global integration point counter
 		int k1,k2,k3;
 		for (int j = 0; j < nGP; j++)
-			for (int i = 0; i < nGP; i++)
+			for (int i = 0; i < nGP; i++, fe.iGP++)
 			{
 				// Local element coordinates and parameter values
 				// of current integration point
@@ -1370,11 +1370,16 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
 					ok = false;
 		}
 
+    // Finalize the element quantities
+    if (ok && !integrand.finalizeElementBou(*A,fe,time))
+      ok = false;
+
 		// Assembly of global system integral
 		if (ok && !glInt.assemble(A->ref(),fe.iel))
 			ok = false;
 		A->destruct();
 
+		firstp += nGP*nGP;
 	}
 
 	return ok;
@@ -1798,63 +1803,18 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol,
 bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
                            const int* npe, char project) const
 {
-	std::cerr << "Fuck it... its the one with the projection in it\n";
+  if (project != 0) {
-	std::cerr << "ASMu3D::evalSolution(...) is not properly implemented yet :(" << std::endl;
+    std::cerr << "ASMu3D::evalSolution() projection schemes not implemented yet\n";
-	return true;
+    return false;
-	// exit(776658);
+  }
-#if 0
-	// Project the secondary solution onto the spline basis
-	Go::SplineVolume* v = nullptr;
-	if (project == 'A')
-		v = this->projectSolutionLocalApprox(integrand);
-	else if (project == 'L')
-		v = this->projectSolutionLocal(integrand);
-	else if (project == 'W')
-		v = this->projectSolutionLeastSquare(integrand);
-	else if (project == 'D' || !npe)
-		v = this->projectSolution(integrand);
 
-	if (npe)
+  // Compute parameter values of the result sampling points
-	{
+  std::array<RealArray,3> gpar;
-		// Compute parameter values of the result sampling points
+  for (int dir = 0; dir < 3; dir++)
-		std::array<RealArray,3> gpar;
+    if (!this->getGridParameters(gpar[dir],dir,npe[dir]-1))
-		if (this->getGridParameters(gpar[0],0,npe[0]-1) &&
+      return false;
-	this->getGridParameters(gpar[1],1,npe[1]-1) &&
+  return this->evalSolution(sField,integrand,gpar.data());
-	this->getGridParameters(gpar[2],2,npe[2]-1))
-		{
-			if (!project)
-				// Evaluate the secondary solution directly at all sampling points
-				return this->evalSolution(sField,integrand,gpar.data());
-			else if (v)
-			{
-				// Evaluate the projected field at the result sampling points
-				const Vector& svec = sField; // using utl::matrix cast operator
-				sField.resize(v->dimension(),
-				              gpar[0].size()*gpar[1].size()*gpar[2].size());
-				v->gridEvaluator(gpar[0],gpar[1],gpar[2],const_cast<Vector&>(svec));
-				delete v;
-				return true;
-			}
-		}
-		else if (v)
-			delete v;
-	}
-	else if (v)
-	{
-		// Extract control point values from the spline object
-		sField.resize(v->dimension(),
-				          v->numCoefs(0)*v->numCoefs(1)*v->numCoefs(2));
-		sField.fill(&(*v->coefs_begin()));
-		delete v;
-		return true;
-	}
 
-	std::cerr <<" *** ASMu3D::evalSolution: Failure!";
-	if (project) std::cerr <<" project="<< project;
-	std::cerr << std::endl;
-	return false;
-#endif
-	return false;
 }