singular dbc corrected

AMDiS/CMakeLists.txt

@@ -497,7 +497,7 @@ if(ENABLE_EXTENSIONS)
     # 	      ${EXTENSIONS_DIR}/base_problems/NavierStokes_Chorin.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesCahnHilliard.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesPhase_TaylorHood.cc
-		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc
+# 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood_RB.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase_RB.cc
@@ -708,6 +708,69 @@ list(APPEND deb_add_dirs "share/amdis/doc")
 
 list(REMOVE_DUPLICATES deb_add_dirs)
 
+# TESTING
+
+ENABLE_TESTING()
+INCLUDE(CTest)
+
+set(TEST_FILE amdis_test.cc)
+file(WRITE ${TEST_FILE} "#include <AMDiS.h>
+int main (int argc, char* argv[]) {
+  AMDiS::init(argc, argv);
+  AMDiS::finalize();
+  return 0 ;
+}")
+
+
+if(WIN32)
+  set(Boost_USE_STATIC_LIBS ON)
+  find_package(Boost 1.42 REQUIRED system iostreams filesystem program_options date_time zlib bzip2)
+else()
+  set(Boost_USE_STATIC_LIBS OFF)
+  find_package(Boost 1.42 REQUIRED system iostreams filesystem program_options date_time)
+endif()
+if(Boost_FOUND)
+  list(APPEND AMDIS_LIBRARIES ${Boost_LIBRARIES})
+  list(APPEND AMDIS_LIBRARY_DIRS ${Boost_LIBRARY_DIRS})
+  list(APPEND AMDIS_INCLUDE_DIRS ${Boost_INCLUDE_DIR})
+else()
+  message(ERROR "Boost libraries not found")
+endif(Boost_FOUND)
+find_library(UMFPACK_LIB umfpack
+	HINTS ENV LIBRARY_PATH
+	DOC "The UMFPACK library")
+if(UMFPACK_LIB)
+	list(APPEND AMDIS_LIBRARIES ${UMFPACK_LIB})
+else()
+	message(FATAL_ERROR "Could not find the UMFPACK library")
+endif()
+
+add_executable("amdis_test" ${TEST_FILE})
+target_link_libraries("amdis_test" amdis reinit compositeFEM muparser ${AMDIS_LIBRARIES})
+
+add_test(NAME test.amdis
+  WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY}
+  COMMAND ${CMAKE_COMMAND} --build . --target amdis_test --config $<CONFIGURATION>
+)
+
+# add_test(NAME test.build.amdis
+#   WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY}
+#   COMMAND ${CMAKE_COMMAND} --build . --target amdis --config $<CONFIGURATION>
+# )
+# add_test(NAME test.build.reinit
+#   WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY}
+#   COMMAND ${CMAKE_COMMAND} --build . --target reinit --config $<CONFIGURATION>
+# )
+# add_test(NAME test.build.compositeFEM
+#   WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY}
+#   COMMAND ${CMAKE_COMMAND} --build . --target compositeFEM --config $<CONFIGURATION>
+# )
+# add_test(NAME test.build.muparser
+#   WORKING_DIRECTORY ${CMAKE_BINARY_DIRECTORY}
+#   COMMAND ${CMAKE_COMMAND} --build . --target muparser --config $<CONFIGURATION>
+# )
+
+
 # PACKAGES
 # ========================================================
 set(CPACK_PACKAGE_NAME "AMDIS")

demo/CMakeLists.txt

@@ -71,6 +71,12 @@ add_executable("vecheat" ${vecheat})
 target_link_libraries("vecheat" ${BASIS_LIBS})
 install(TARGETS vecheat RUNTIME DESTINATION bin)
 
+
+set(test src/test.cc)
+add_executable("test" ${test})
+target_link_libraries("test" ${BASIS_LIBS})
+install(TARGETS test RUNTIME DESTINATION bin)
+
 # add_executable("matrix_test" src/matrix_test.cc)
 # target_link_libraries("matrix_test" ${BASIS_LIBS})
 # install(TARGETS vecheat RUNTIME DESTINATION bin)

demo/init/ellipt.dat.3d

@@ -15,17 +15,17 @@ ellipt->solver->info:           2
 ellipt->solver->left precon:    diag
 ellipt->solver->right precon:   no
 
-ellipt->estimator[0]:              residual
+%ellipt->estimator[0]:              residual
 ellipt->estimator[0]->error norm:  H1_NORM   % 1: H1_NORM, 2: L2_NORM
 ellipt->estimator[0]->C1:          0.01 % constant of jump residual
 
-ellipt->marker[0]->strategy:       2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+ellipt->marker[0]->strategy:      0 % 2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
 ellipt->marker[0]->MSGamma:        0.5
 
 ellipt->adapt[0]->tolerance:          1e-4
 ellipt->adapt[0]->refine bisections:  2
 
-ellipt->adapt->max iteration:         10
+ellipt->adapt->max iteration:       0 %  10
 
 ellipt->output->filename:        output/ellipt.3d
 ellipt->output->ParaView format: 1

demo/src/ellipt.cc

@@ -62,17 +62,17 @@ int main(int argc, char* argv[])
   
   // ===== create matrix operator =====
   Operator matrixOperator(ellipt.getFeSpace());
-  matrixOperator.addTerm(new Simple_SOT);
+  matrixOperator.addTerm(new Simple_ZOT);
   ellipt.addMatrixOperator(matrixOperator, 0, 0);
 
 
   // ===== create rhs operator =====
-  Operator rhsOperator(ellipt.getFeSpace());
-  rhsOperator.addTerm(new CoordsAtQP_ZOT(new F));
-  ellipt.addVectorOperator(rhsOperator, 0);
+//   Operator rhsOperator(ellipt.getFeSpace());
+//   rhsOperator.addTerm(new CoordsAtQP_ZOT(new F));
+//   ellipt.addVectorOperator(rhsOperator, 0);
   
   // ===== add boundary conditions =====
-  ellipt.addDirichletBC(1, 0, 0, new G);
+  ellipt.addDirichletBC(1, 0, 0, new AMDiS::Const<double, WorldVector<double> >(1.0));
 
 
   // ===== start adaption loop =====

extensions/ExtendedProblemStat.h

@@ -68,11 +68,15 @@ public:
   
   ~ExtendedProblemStat()
   {
-    for (int i = 0; i < nComponents; ++i)
+    for (int i = 0; i < nComponents; ++i) {
       if (exactSolutions[i]) {
 	delete exactSolutions[i];
 	exactSolutions[i] = NULL;
       }
+    }
+    
+    for (size_t i = 0; i < DirichletBcDataList.size(); i++)
+      delete DirichletBcDataList[i];
     
   }
 
@@ -80,7 +84,7 @@ public:
   void initialize(Flag initFlag,
 		  ProblemStatSeq *adoptProblem = NULL,
 		  Flag adoptFlag = INIT_NOTHING)
-  {
+  {    
     ProblemStat_::initialize(initFlag, adoptProblem, adoptFlag);
     for (int i = 0; i < nComponents; ++i)
       exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution");
@@ -144,7 +148,7 @@ public:
     // apply dirichlet boundary conditions
     for (size_t k = 0; k < singularDirichletBC.size(); k++)
       applyDirichletBC(singularDirichletBC[k], asmMatrix, asmVector);
-
+    
     // update solverMatrix
     if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
       solverMatrix.setMatrix(*getSystemMatrix());
@@ -230,7 +234,7 @@ public:
   {
     WARNING("Depreciated! Use addDirichletBC instead!!!\n");
     DofIndex* dofIndex = new DofIndex(idx);
-    addDirichletBC(*dofIndex, row, col, values);
+    addDirichletBC(dofIndex, row, col, &values);
   }
 
   
@@ -333,7 +337,9 @@ public:
     
 protected:
   // traverse matrix rows and set unity row where dirichlet condition shall be applied.
-  void applyDirichletBC(size_t row_, size_t col_, DegreeOfFreedom idx, double value, bool asmMatrix = true, bool asmVector = true)
+  void applyDirichletBC(size_t row_, size_t col_, 
+			DegreeOfFreedom idx, double value, 
+			bool asmMatrix = true, bool asmVector = true)
   {
     using namespace mtl;
     typedef DOFMatrix::base_matrix_type Matrix;

extensions/ProblemStatMassConserve.h

@@ -125,7 +125,7 @@ struct ProblemStatMassConserve : public ExtendedProblemStat
     prob2->setMeshes(meshes2);    
     prob2->initialize(INIT_ALL - INIT_MESH);
     
-    for (size_t i = 1; i < meshes2.size(); i++) {
+    for (size_t i = 0; i < meshes2.size(); i++) {
       int globalRefinements = 0;
 
       // If AMDiS is compiled for parallel computations, the global refinements are
@@ -139,12 +139,10 @@ struct ProblemStatMassConserve : public ExtendedProblemStat
       bool initMesh = true;
 
       // Initialize the meshes if there is no serialization file.
-      if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized()))
+      if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized())) {
 	meshes2[i]->initialize();
-
-      // do global refinements
-      if (initMesh && meshes2[i])
 	prob2->getRefinementManager()->globalRefine(meshes2[i], globalRefinements);
+      }
     }
     
     

extensions/SingularDirichletBC2.cc

@@ -43,7 +43,7 @@ void getIndicesAux<WorldVector<double> >(
   
 // AbstractFunction (signedDistFct)
 template<>
-void getIndicesAux(
+void getIndicesAux<AbstractFunction<double, WorldVector<double> > >(
 	AbstractFunction<double, WorldVector<double> >& pos, 
 	BoundaryType boundary_nr,
 	const FiniteElemSpace *feSpace, 
@@ -67,7 +67,7 @@ void getIndicesAux<DOFVector<double> >(
 
 // Meshindicator
 template<>
-void getIndicesAux(
+void getIndicesAux<AbstractFunction<bool, WorldVector<double> > >(
 	AbstractFunction<bool, WorldVector<double> >& pos, 
 	BoundaryType boundary_nr,
 	const FiniteElemSpace *feSpace, 
@@ -85,7 +85,8 @@ void getValuesAux<double>(
 	const double &val, 
 	const FiniteElemSpace *feSpace, 
 	std::map<DegreeOfFreedom, double> &indicesValues,
-	boost::false_type is_not_abstract_fct)
+	boost::false_type is_not_abstract_fct,
+	int depth)
 {
   std::map<DegreeOfFreedom, double>::iterator it;
   for (it = indicesValues.begin(); it != indicesValues.end(); it++)
@@ -97,13 +98,14 @@ void getValuesAux<AbstractFunction<double, WorldVector<double> > >(
 	const AbstractFunction<double, WorldVector<double> >& val, 
 	const FiniteElemSpace *feSpace, 
 	std::map<DegreeOfFreedom, double> &indicesValues,
-	boost::true_type is_abstract_fct)
+	boost::true_type is_abstract_fct,
+	int depth)
 {
-  WorldVector<double> x;
+  DOFVector<WorldVector<double> > x(feSpace, "coords");
+  feSpace->getMesh()->getDofIndexCoords(x);
   std::map<DegreeOfFreedom, double>::iterator it;
   for (it = indicesValues.begin(); it != indicesValues.end(); it++)
-    if (feSpace->getMesh()->getDofIndexCoords(it->first, feSpace, x))
-      it->second = val(x);
+    it->second = val(x[it->first]);
 }
 
 // DOFVector<double>
@@ -112,7 +114,8 @@ void getValuesAux<DOFVector<double> >(
 	const DOFVector<double>& val, 
 	const FiniteElemSpace *feSpace, 
 	std::map<DegreeOfFreedom, double> &indicesValues,
-	boost::false_type is_not_abstract_fct)
+	boost::false_type is_not_abstract_fct,
+	int depth)
 {
   std::map<DegreeOfFreedom, double>::iterator it;
   for (it = indicesValues.begin(); it != indicesValues.end(); it++)

extensions/SingularDirichletBC2.h

@@ -86,6 +86,23 @@ namespace details {
 		  const FiniteElemSpace* 
 		  feSpace, std::map<DegreeOfFreedom, double> &indices, 
 		  boost::true_type is_abstract_fct);
+  
+  // template specialization for AbstractFunctions
+  template<>
+  void getIndicesAux<AbstractFunction<double, WorldVector<double> > >(
+		  AbstractFunction<double, WorldVector<double> > &pos, 
+		  BoundaryType b, 
+		  const FiniteElemSpace* 
+		  feSpace, std::map<DegreeOfFreedom, double> &indices, 
+		  boost::true_type is_abstract_fct);
+  template<>
+  void getIndicesAux<AbstractFunction<bool, WorldVector<double> > >(
+		  AbstractFunction<bool, WorldVector<double> > &pos, 
+		  BoundaryType b, 
+		  const FiniteElemSpace* 
+		  feSpace, std::map<DegreeOfFreedom, double> &indices, 
+		  boost::true_type is_abstract_fct);
+  
   template<typename PosType>
   void getIndicesAux(PosType &pos, 
 		  BoundaryType b, 
@@ -114,20 +131,34 @@ namespace details {
   void getValuesAux(const ValueType &val, 
 		 const FiniteElemSpace* feSpace, 
 		 std::map<DegreeOfFreedom, double> &indices, 
-		 boost::true_type is_abstract_fct);
+		 boost::true_type is_abstract_fct,
+		 int depth = 0)
+  {
+    TEST_EXIT(depth == 0)("Recursion bricht nicht ab!\n");
+    details::getValuesAux(dynamic_cast<const AbstractFunction<double, WorldVector<double> >&>(val), feSpace, indices, is_abstract_fct, depth+1);
+  }
+  // template specialization for AbstractFunctions
+  template<>
+  void getValuesAux<AbstractFunction<double, WorldVector<double> > >(
+		  const AbstractFunction<double, WorldVector<double> >& val, 
+		  const FiniteElemSpace *feSpace, 
+		  std::map<DegreeOfFreedom, double> &indicesValues,
+		  boost::true_type is_abstract_fct,
+		  int depth);
+    
   template<typename ValueType>
   void getValuesAux(const ValueType &val, 
 		 const FiniteElemSpace* feSpace, 
 		 std::map<DegreeOfFreedom, double> &indices, 
-		 boost::false_type is_not_abstract_fct);
+		 boost::false_type is_not_abstract_fct,
+		 int depth = 0);
 
   template<typename ValueType>
   void getValues(const ValueType &val, 
 		 const FiniteElemSpace* feSpace, 
 		 std::map<DegreeOfFreedom, double> &indices)
   {
-    details::getValuesAux(val, feSpace, indices, 
-	typename boost::is_base_of<AbstractFunction<double, WorldVector<double> >, ValueType>::type());
+    details::getValuesAux(val, feSpace, indices, typename boost::is_base_of<AbstractFunction<double, WorldVector<double> >, ValueType>::type());
   }
 
   /**
@@ -300,7 +331,7 @@ struct DirichletBcData : DirichletBcDataBase {
   DirichletBcData(int i_, int j_, BoundaryType nr_, PosType &meshIndicator_, const ValueType &val_)
   : DirichletBcDataBase(i_, j_), boundary_nr(nr_), pos(meshIndicator_), val(val_) {}
     
-  virtual void addToList(const FiniteElemSpace *feSpace, std::vector<SingularDirichletBC> &list)
+  void addToList(const FiniteElemSpace *feSpace, std::vector<SingularDirichletBC> &list) override
   {
     std::map<DegreeOfFreedom, double> indicesValues;
     

extensions/base_problems/CahnHilliard.cc

@@ -24,7 +24,10 @@
 
 using namespace AMDiS;
 
-CahnHilliard::CahnHilliard(const std::string &name_) :
+namespace details {
+
+template<typename P>
+CahnHilliard<P>::CahnHilliard(const std::string &name_) :
   super(name_),
   useMobility(false),
   useReinit(false),
@@ -45,7 +48,7 @@ CahnHilliard::CahnHilliard(const std::string &name_) :
   // type of double well: 0= [0,1], 1= [-1,1]
   Parameters::get(name_ + "->double-well type", doubleWell); 
   
-  Parameters::get(name + "->use reinit", useReinit);
+  Parameters::get(name_ + "->use reinit", useReinit);
 
   // transformation of the parameters
   minusEps = -eps;
@@ -56,106 +59,106 @@ CahnHilliard::CahnHilliard(const std::string &name_) :
 }
 
 
-void CahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo) 
-{ FUNCNAME("CahnHilliard::solveInitialProblem()");
-
-  Flag initFlag = initDataFromFile(adaptInfo);
+template<typename P>
+void CahnHilliard<P>::solveInitialProblem(AdaptInfo *adaptInfo) 
+{
+  Flag initFlag = self::initDataFromFile(adaptInfo);
 
   if (!initFlag.isSet(DATA_ADOPTED)) {
     int initialInterface = 0;
-    Initfile::get(name + "->initial interface", initialInterface);
+    Initfile::get(self::name + "->initial interface", initialInterface);
     double initialEps = eps;
-    Initfile::get(name + "->initial epsilon", initialEps);
+    Initfile::get(self::name + "->initial epsilon", initialEps);
 
     if (initialInterface == 0) {
       /// horizontale Linie
       double a= 0.0, dir= -1.0;
-      Initfile::get(name + "->line->pos", a);
-      Initfile::get(name + "->line->direction", dir);
-      prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
+      Initfile::get(self::name + "->line->pos", a);
+      Initfile::get(self::name + "->line->direction", dir);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
     }
     else if (initialInterface == 1) {
       /// schraege Linie
       double theta = m_pi/4.0;
-      prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0));
-      transformDOFInterpolation(prob->getSolution()->getDOFVector(0),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min<double>);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0));
+      transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min<double>);
     }
     else if (initialInterface == 2) {
       /// Ellipse
       double a= 1.0, b= 1.0;
-      Initfile::get(name + "->ellipse->a", a);
-      Initfile::get(name + "->ellipse->b", b);
-      prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b));
+      Initfile::get(self::name + "->ellipse->a", a);
+      Initfile::get(self::name + "->ellipse->b", b);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b));
     }
     else if (initialInterface == 3) {
       /// zwei horizontale Linien
       double a= 0.75, b= 0.375;
-      Initfile::get(name + "->lines->pos1", a);
-      Initfile::get(name + "->lines->pos2", b);
-      prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0));
-      transformDOFInterpolation(prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max<double>);
+      Initfile::get(self::name + "->lines->pos1", a);
+      Initfile::get(self::name + "->lines->pos2", b);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0));
+      transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max<double>);
     }
     else if (initialInterface == 4) {
       /// Kreis
       double radius= 1.0;
-      Initfile::get(name + "->kreis->radius", radius);
-      prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius));
+      Initfile::get(self::name + "->kreis->radius", radius);
+      self::prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius));
     } else if (initialInterface == 5) {
       /// Rechteck
       double width = 0.5;
       double height = 0.3;
       WorldVector<double> center; center.set(0.5);
-      Initfile::get(name + "->rectangle->width", width);
-      Initfile::get(name + "->rectangle->height", height);
-      Initfile::get(name + "->rectangle->center", center);
-      prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center));
+      Initfile::get(self::name + "->rectangle->width", width);
+      Initfile::get(self::name + "->rectangle->height", height);
+      Initfile::get(self::name + "->rectangle->center", center);
+      self::prob->getSolution()->getDOFVector(0)->inte