transformDOF extended

AMDiS/src/TransformDOF.h

@@ -547,12 +547,142 @@ inline void transformDOF(DOFVector<T1> *vec1, DOFVector<T2> *vec2, T3 val, DOFVe
     }
 }
 
-// ===========================================================================================
-
 template<typename T1, typename T2, typename T3, typename T4>
 inline void transformDOF(DOFVector<T1> &vec1, DOFVector<T2> &vec2, T3 val, DOFVector<T4> &result, TertiaryAbstractFunction<T4, T1, T2, T3> &tertiary_op)
 { transformDOF(&vec1, &vec2, val, &result, &tertiary_op); }
-	
+
+// result = tertiary_op(vec1, vec2, vec3)
+template<typename T1, typename T2, typename T3, typename T4>
+inline void transformDOF(DOFVector<T1> *vec1, T2 val, DOFVector<T2> *vec3, DOFVector<T4> *result, TertiaryAbstractFunction<T4, T1, T2, T3> *tertiary_op)
+{
+
+    DOFVector<T4>* res;
+    bool useResult = true;
+    if(static_cast<void*>(vec1)==static_cast<void*>(result)
+      || static_cast<void*>(vec3) == static_cast<void*>(result)) {
+      res= new DOFVector<T4>(result->getFeSpace(), result->getName());
+      useResult = false;
+    } else
+      res= result;
+
+    const FiniteElemSpace *vec1FeSpace = vec1->getFeSpace();
+    const FiniteElemSpace *vec3FeSpace = vec3->getFeSpace();
+    const FiniteElemSpace *resFeSpace = res->getFeSpace();
+
+    const BasisFunction *vec1BasisFcts = vec1FeSpace->getBasisFcts();
+    const BasisFunction *vec3BasisFcts = vec3FeSpace->getBasisFcts();
+    const BasisFunction *resBasisFcts = resFeSpace->getBasisFcts();
+
+    int nVec1BasisFcts = vec1BasisFcts->getNumber();
+    int nVec3BasisFcts = vec3BasisFcts->getNumber();
+    int nResBasisFcts = resBasisFcts->getNumber();
+
+    std::vector<DegreeOfFreedom> resLocalIndices(nResBasisFcts);
+    ElementVector vec1LocalCoeffs(nVec1BasisFcts);
+    ElementVector vec3LocalCoeffs(nVec3BasisFcts);
+
+    DimVec<double> *coords1 = NULL, *coords2 = NULL;
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(resFeSpace->getMesh(), -1,
+            Mesh::CALL_LEAF_EL |
+            Mesh::FILL_COORDS);
+
+    while (elInfo) {
+      Element *el = elInfo->getElement();
+
+      resBasisFcts->getLocalIndices(el, resFeSpace->getAdmin(), resLocalIndices);
+      vec1->getLocalVector(el, vec1LocalCoeffs);
+      vec3->getLocalVector(el, vec3LocalCoeffs);
+
+      for (int i = 0; i < nResBasisFcts; i++) {
+        coords1 = vec1BasisFcts->getCoords(i);
+        coords2 = vec3BasisFcts->getCoords(i);
+        (*res)[resLocalIndices[i]] = (*tertiary_op)(
+          vec1BasisFcts->evalUh(*coords1, vec1LocalCoeffs),
+          val,
+          vec3BasisFcts->evalUh(*coords2, vec3LocalCoeffs)
+	  );
+      }
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    if (!useResult) {
+      result->copy(*res);
+      delete res;
+    }
+}
+
+template<typename T1, typename T2, typename T3, typename T4>
+inline void transformDOF(DOFVector<T1> &vec1, T2 val, DOFVector<T3> &vec3, DOFVector<T4> &result, TertiaryAbstractFunction<T4, T1, T2, T3> &tertiary_op)
+{ transformDOF(&vec1, val, &vec3, &result, &tertiary_op); }
+
+// result = tertiary_op(vec1, vec2, vec3)
+template<typename T1, typename T2, typename T3, typename T4>
+inline void transformDOF(T1 val, DOFVector<T2> *vec2, DOFVector<T2> *vec3, DOFVector<T4> *result, TertiaryAbstractFunction<T4, T1, T2, T3> *tertiary_op)
+{
+
+    DOFVector<T4>* res;
+    bool useResult = true;
+    if(static_cast<void*>(vec2)==static_cast<void*>(result)
+      || static_cast<void*>(vec3) == static_cast<void*>(result)) {
+      res= new DOFVector<T4>(result->getFeSpace(), result->getName());
+      useResult = false;
+    } else
+      res= result;
+
+    const FiniteElemSpace *vec2FeSpace = vec2->getFeSpace();
+    const FiniteElemSpace *vec3FeSpace = vec3->getFeSpace();
+    const FiniteElemSpace *resFeSpace = res->getFeSpace();
+
+    const BasisFunction *vec2BasisFcts = vec2FeSpace->getBasisFcts();
+    const BasisFunction *vec3BasisFcts = vec3FeSpace->getBasisFcts();
+    const BasisFunction *resBasisFcts = resFeSpace->getBasisFcts();
+
+    int nVec2BasisFcts = vec2BasisFcts->getNumber();
+    int nVec3BasisFcts = vec3BasisFcts->getNumber();
+    int nResBasisFcts = resBasisFcts->getNumber();
+
+    std::vector<DegreeOfFreedom> resLocalIndices(nResBasisFcts);
+    ElementVector vec2LocalCoeffs(nVec2BasisFcts);
+    ElementVector vec3LocalCoeffs(nVec3BasisFcts);
+
+    DimVec<double> *coords1 = NULL, *coords2 = NULL;
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(resFeSpace->getMesh(), -1,
+            Mesh::CALL_LEAF_EL |
+            Mesh::FILL_COORDS);
+
+    while (elInfo) {
+      Element *el = elInfo->getElement();
+
+      resBasisFcts->getLocalIndices(el, resFeSpace->getAdmin(), resLocalIndices);
+      vec2->getLocalVector(el, vec2LocalCoeffs);
+      vec3->getLocalVector(el, vec3LocalCoeffs);
+
+      for (int i = 0; i < nResBasisFcts; i++) {
+        coords1 = vec2BasisFcts->getCoords(i);
+        coords2 = vec3BasisFcts->getCoords(i);
+        (*res)[resLocalIndices[i]] = (*tertiary_op)(
+          val,
+          vec2BasisFcts->evalUh(*coords1, vec2LocalCoeffs),
+          vec3BasisFcts->evalUh(*coords2, vec3LocalCoeffs)
+	  );
+      }
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    if (!useResult) {
+      result->copy(*res);
+      delete res;
+    }
+}
+
+template<typename T1, typename T2, typename T3, typename T4>
+inline void transformDOF(T1 val, DOFVector<T2> &vec2, DOFVector<T3> &vec3, DOFVector<T4> &result, TertiaryAbstractFunction<T4, T1, T2, T3> &tertiary_op)
+{ transformDOF(val, &vec2, &vec3, &result, &tertiary_op); }
+
+// ===========================================================================================
+
 // return binary_op(vec, interpol(fct))
 template<typename T>
 inline void transformDOFInterpolation(