some new baseproblems

extensions/ExtendedProblemStat.h

@@ -11,9 +11,9 @@
 
 using namespace AMDiS;
 
-const Flag INIT_EXACT_SOLUTION = 0X2000L;
-const Flag UPDATE_PERIODIC_BC = 0X4000L;
-const Flag UPDATE_DIRICHLET_BC = 0X8000L;
+const Flag INIT_EXACT_SOLUTION = 0X10000L;
+const Flag UPDATE_PERIODIC_BC  = 0X20000L;
+const Flag UPDATE_DIRICHLET_BC = 0X40000L;
 
 #if defined NONLIN_PROBLEM
 typedef ProblemNonLin ProblemStat_;
@@ -128,7 +128,10 @@ public:
     if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
       solverMatrix.setMatrix(*getSystemMatrix());
     }
-//     writeMatrix("matrix.mat");
+#endif
+
+#ifndef NDEBUG
+     writeMatrix(name + ".mtx");
 #endif
   }
 

extensions/Helpers.h

@@ -205,7 +205,7 @@ namespace Helpers {
   
   /// calculate the dimension of a mesh, by mesh traversal.
   void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner);
-
+  
   /// read DOFVector from AMDiS .dat-files
   void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size);
 

extensions/POperators.cc

@@ -113,27 +113,27 @@ WorldVector<double> Phase_FOT::f(const int iq) const
 Phase_ZOT::Phase_ZOT(DOFVectorBase<double>* phaseDV_, double fac_) : ZeroOrderTerm(6),
 	phaseDV(phaseDV_), fac(fac_)
 {	
-	TEST_EXIT(phaseDV_)("No vector Phase!\n");
+  TEST_EXIT(phaseDV_)("No vector Phase!\n");
 	
-	auxFeSpaces.insert(phaseDV_->getFeSpace());
+  auxFeSpaces.insert(phaseDV_->getFeSpace());
 }
 void Phase_ZOT::initElement(const ElInfo* elInfo, 
 			SubAssembler* subAssembler,
 			Quadrature *quad) 
 {
-	getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
+  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
 }
 void Phase_ZOT::initElement(const ElInfo* largeElInfo,
 			const ElInfo* smallElInfo,
 			SubAssembler* subAssembler,
 			Quadrature *quad)
 {
-	getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
+  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
 }
 void Phase_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
 { 
-	for (int iq = 0; iq < nPoints; iq++)
-		C[iq] += f(iq) * phase[iq] * fac;
+  for (int iq = 0; iq < nPoints; iq++)
+    C[iq] += f(iq) * phase[iq] * fac;
 }
 void Phase_ZOT::eval(int nPoints,
 	const mtl::dense_vector<double>& uhAtQP,
@@ -142,13 +142,13 @@ void Phase_ZOT::eval(int nPoints,
 	mtl::dense_vector<double>& result,
 	double opFactor)
 {
-	double factor = opFactor * fac;
-	for (int iq = 0; iq < nPoints; iq++)
-		result[iq] += factor * f(iq) * phase[iq] * uhAtQP[iq];
+  double factor = opFactor * fac;
+  for (int iq = 0; iq < nPoints; iq++)
+    result[iq] += factor * f(iq) * phase[iq] * uhAtQP[iq];
 }
 double Phase_ZOT::f(const int iq) const
 {
-	return 1.0;
+  return 1.0;
 }
 
 /* ----------------------------------------------------------- */
@@ -2586,3 +2586,42 @@ double Viscosity2_SOT::f(const int iq) const
     double viscosity= v1*phase+v2*(1.0-phase);
     return viscosity;
 }
+
+
+
+  // =========== VecGrad2_FOT ===========
+
+  void VecGrad2_FOT::initElement(const ElInfo* elInfo, 
+				SubAssembler* subAssembler,
+				Quadrature *quad) 
+  {
+    getVectorAtQPs(vec0, elInfo, subAssembler, quad, vec0AtQPs);
+    getGradientsAtQPs(vec1, elInfo, subAssembler, quad, grad1AtQPs);
+    getGradientsAtQPs(vec2, elInfo, subAssembler, quad, grad2AtQPs);
+  }
+
+  void VecGrad2_FOT::getLb(const ElInfo *elInfo,
+			   vector<mtl::dense_vector<double> >& Lb) const 
+  {
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+    const int nPoints = static_cast<int>(Lb.size());
+
+    for (int iq = 0; iq < nPoints; iq++)
+      lb(grdLambda, (*vecFct)(vec0AtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
+  }
+  
+  void VecGrad2_FOT::eval(int nPoints,
+			  const mtl::dense_vector<double>& uhAtQP,
+			  const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+			  const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+			  mtl::dense_vector<double>& result,
+			  double factor) 
+  {
+    if (num_rows(grdUhAtQP) > 0) {
+      for (int iq = 0; iq < nPoints; iq++) {
+	WorldVector<double> b = (*vecFct)(vec0AtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
+	result[iq] += b * grdUhAtQP[iq] * factor;
+      }
+    }
+  }
+

extensions/POperators_FOT.h

@@ -546,6 +546,70 @@ mtl::dense_vector<double> ch;
 mtl::dense_vector<WorldVector<double> > grdMu;
 double delta;
 
+
 };
 
+/* -------------------------------------------------------------- */
+
+
+
+  /**
+   * \ingroup Assembler
+   *
+   * \brief
+   * First order term: \f$ b(v(\vec{x}), \nabla w_1(\vec{x}), \nabla w_2(\vec{x})) \cdot \nabla u(\vec{x})\f$.
+   */
+  class VecGrad2_FOT : public FirstOrderTerm
+  {
+  public:
+    /// Constructor
+    VecGrad2_FOT(DOFVectorBase<double> *dv, DOFVectorBase<double> *dw1, DOFVectorBase<double> *dw2,
+		TertiaryAbstractFunction<WorldVector<double>, double, WorldVector<double>, WorldVector<double> > *vecFct_)
+      : FirstOrderTerm(vecFct_->getDegree()), vec0(dv), vec1(dw1), vec2(dw2), vecFct(vecFct_) 
+    
+    {   
+      auxFeSpaces.insert(dv->getFeSpace()); 
+      auxFeSpaces.insert(dw1->getFeSpace()); 
+      auxFeSpaces.insert(dw2->getFeSpace());
+    }
+
+    /// Implementation of \ref OperatorTerm::initElement().
+    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+		     Quadrature *quad = NULL);
+    
+    /// Implements FirstOrderTerm::getLb().
+    void getLb(const ElInfo *elInfo,
+	       vector<mtl::dense_vector<double> >& Lb) const; 
+    
+    /// Implements FirstOrderTerm::eval().
+    void eval(int nPoints,
+	      const mtl::dense_vector<double>& uhAtQP,
+	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+	      mtl::dense_vector<double>& result,
+	      double factor);
+    
+  protected:
+    /// DOFVector to be evaluated at quadrature points.
+    DOFVectorBase<double>* vec0;
+    
+    /// DOFVector to be evaluated at quadrature points.
+    DOFVectorBase<double>* vec1;
+
+    /// DOFVector to be evaluated at quadrature points.
+    DOFVectorBase<double>* vec2;
+    
+    mtl::dense_vector<double> vec0AtQPs;
+    
+    /// Gradient v at quadrature points.
+    mtl::dense_vector<WorldVector<double> > grad1AtQPs;
+
+    /// Gradient v at quadrature points.
+    mtl::dense_vector<WorldVector<double> > grad2AtQPs;
+
+    /// Function for b.
+    TertiaryAbstractFunction<WorldVector<double>, double, WorldVector<double>, WorldVector<double> > *vecFct;
+  };
+
+
 #endif // P_OPERATORS_FOT_H