ln elast dist

.gitignore

@@ -4,7 +4,13 @@
 # in that subdirectory instead.
 #
 
+*.out.*
+*.tmp
+*.cbp
+*.m
+.gitignore
 CMake.vars
+CTestCheckpoint.txt
 *.h2d
 .*
 *.o
@@ -18,7 +24,6 @@ tutorial.*
 *.png
 *.bmp
 *.user
-*make*
 *.bak
 
 # top-level files
@@ -65,9 +70,6 @@ Debug/
 # kdevelop project files
 *.kdev*
 
-# Python stuff
-!python/_hermes_common.so
-
 # Documentation generated
 *.html
 *.eps
@@ -80,4 +82,21 @@ Debug/
 *.css
 *.tag
 *.sty
-*.dvi
+*.dvi
+*.vtk
+
+
+# Misc
+*.exe
+*.dat
+*.lastbuildstate
+*.cfg
+*.bak
+*.options
+*.pdr.*
+*.jit
+*My Inspector*
+*.exp
+*.ipgset
+*.psess
+*.vsp

1d/layer-boundary/definitions.cpp

@@ -15,7 +15,7 @@ double CustomExactFunction::dduhat_dxx(double x)
   return -K*K * (exp(K*x) + exp(-K*x)) / (exp(K) + exp(-K));
 }
 
-CustomExactSolution::CustomExactSolution(const Mesh* mesh, double K) : ExactSolutionScalar<double>(mesh)
+CustomExactSolution::CustomExactSolution(MeshSharedPtr mesh, double K) : ExactSolutionScalar<double>(mesh), K(K)
 {
   cef = new CustomExactFunction(K);
 }
@@ -36,7 +36,7 @@ void CustomExactSolution::derivatives(double x, double y, double& dx, double& dy
   dy = 0;
 }
 
-Ord CustomExactSolution::ord(Ord x, Ord y) const 
+Ord CustomExactSolution::ord(double x, double y) const 
 {
   return Ord(20);
 }

1d/layer-boundary/definitions.h

@@ -25,17 +25,20 @@ class CustomExactFunction
 class CustomExactSolution : public ExactSolutionScalar<double>
 {
 public:
-  CustomExactSolution(const Mesh* mesh, double K);
+  CustomExactSolution(MeshSharedPtr mesh, double K);
 
   virtual double value (double x, double y) const;
 
   virtual void derivatives(double x, double y, double& dx, double& dy) const;
 
-  virtual Ord ord(Ord x, Ord y) const;
+  virtual Ord ord (double x, double y) const;
+
+  MeshFunction<double>* clone() const { return new CustomExactSolution(this->mesh, this->K); }
 
   ~CustomExactSolution();
 
   CustomExactFunction* cef;
+  double K;
 };
 
 /* Custom function */

1d/layer-boundary/main.cpp

@@ -1,5 +1,5 @@
-#define HERMES_REPORT_ALL
-#define HERMES_REPORT_FILE "application.log"
+
+
 #include "definitions.h"
 
 using namespace RefinementSelectors;
@@ -24,18 +24,13 @@ const int P_INIT = 1;
 const int INIT_REF_NUM = 0;                       
 // Number of initial mesh refinements towards the boundary.
 const int INIT_REF_NUM_BDY = 5;                   
-// This is a quantitative parameter of the adapt(...) function and
-// it has different meanings for various adaptive strategies (see below).
-const double THRESHOLD = 0.3;                     
-// Adaptive strategy:
-// STRATEGY = 0 ... refine elements until sqrt(THRESHOLD) times total
-//   error is processed. If more elements have similar errors, refine
-//   all to keep the mesh symmetric.
-// STRATEGY = 1 ... refine all elements whose error is larger
-//   than THRESHOLD times maximum element error.
-// STRATEGY = 2 ... refine all elements whose error is larger
-//   than THRESHOLD.
-const int STRATEGY = 0;                           
+/// This is a quantitative parameter of the adapt(...) function and
+// it has different meanings for various adaptive strategies.
+const double THRESHOLD = 0.5;
+// Error calculation & adaptivity.
+DefaultErrorCalculator<double, HERMES_H1_NORM> errorCalculator(RelativeErrorToGlobalNorm, 1);
+// Stopping criterion for an adaptivity step.
+AdaptStoppingCriterionSingleElement<double> stoppingCriterion(THRESHOLD);
 // Predefined list of element refinement candidates. Possible values are
 // H2D_P_ISO, H2D_P_ANISO, H2D_H_ISO, H2D_H_ANISO, H2D_HP_ISO,
 // H2D_HP_ANISO_H, H2D_HP_ANISO_P, H2D_HP_ANISO.
@@ -66,19 +61,19 @@ const double K = 1e2;
 int main(int argc, char* argv[])
 {
   // Load the mesh.
-  Mesh mesh;
+  MeshSharedPtr mesh(new Mesh);
   MeshReaderH1DXML mloader;
-  mloader.load("domain.xml", &mesh);
+  mloader.load("domain.xml", mesh);
 
   // Perform initial mesh refinement.
   // Split elements vertically.
   int refinement_type = 2;                        
-  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements(refinement_type);
-  mesh.refine_towards_boundary("Left", INIT_REF_NUM_BDY);
-  mesh.refine_towards_boundary("Right", INIT_REF_NUM_BDY);
+  for (int i = 0; i < INIT_REF_NUM; i++) mesh->refine_all_elements(refinement_type);
+  mesh->refine_towards_boundary("Left", INIT_REF_NUM_BDY);
+  mesh->refine_towards_boundary("Right", INIT_REF_NUM_BDY);
 
   // Define exact solution.
-  CustomExactSolution exact_sln(&mesh, K);
+  MeshFunctionSharedPtr<double> exact_sln(new CustomExactSolution(mesh, K));
 
   // Define right side vector.
   CustomFunction f(K);
@@ -91,13 +86,13 @@ int main(int argc, char* argv[])
   EssentialBCs<double> bcs(&bc_essential);
 
   // Create an H1 space with default shapeset.
-  H1Space<double> space(&mesh, &bcs, P_INIT);
+  SpaceSharedPtr<double> space(new H1Space<double>(mesh, &bcs, P_INIT));
 
   // Initialize approximate solution.
-  Solution<double> sln;
+  MeshFunctionSharedPtr<double> sln(new Solution<double>());
 
   // Initialize refinement selector.
-  H1ProjBasedSelector<double> selector(CAND_LIST, CONV_EXP, H2DRS_DEFAULT_ORDER);
+  H1ProjBasedSelector<double> selector(CAND_LIST);
 
   // Initialize views.
   Views::ScalarView sview("Solution", new Views::WinGeom(0, 0, 600, 360));
@@ -123,11 +118,11 @@ int main(int argc, char* argv[])
     int order_increase = 1;          
     // FIXME: This should be '2' but that leads to a segfault.
     int refinement_type = 0;         
-    Mesh::ReferenceMeshCreator refMeshCreator(&mesh);
-    Mesh* ref_mesh = refMeshCreator.create_ref_mesh();
+    Mesh::ReferenceMeshCreator refMeshCreator(mesh);
+    MeshSharedPtr ref_mesh = refMeshCreator.create_ref_mesh();
 
-    Space<double>::ReferenceSpaceCreator refSpaceCreator(&space, ref_mesh);
-    Space<double>* ref_space = refSpaceCreator.create_ref_space();
+    Space<double>::ReferenceSpaceCreator refSpaceCreator(space, ref_mesh);
+    SpaceSharedPtr<double> ref_space = refSpaceCreator.create_ref_space();
     int ndof_ref = ref_space->get_num_dofs();
 
     Hermes::Mixins::Loggable::Static::info("---- Adaptivity step %d (%d DOF):", as, ndof_ref);
@@ -141,7 +136,7 @@ int main(int argc, char* argv[])
     NewtonSolver<double> newton(&dp);
     //newton.set_verbose_output(false);
 
-    Solution<double> ref_sln;
+    MeshFunctionSharedPtr<double> ref_sln(new Solution<double>());
     try
     {
       newton.solve();
@@ -152,43 +147,45 @@ int main(int argc, char* argv[])
       throw Hermes::Exceptions::Exception("Newton's iteration failed.");
     };
     // Translate the resulting coefficient vector into the instance of Solution.
-    Solution<double>::vector_to_solution(newton.get_sln_vector(), ref_space, &ref_sln);
+    Solution<double>::vector_to_solution(newton.get_sln_vector(), ref_space, ref_sln);
 
     cpu_time.tick();
     Hermes::Mixins::Loggable::Static::info("Solution: %g s", cpu_time.last());
 
     // Project the fine mesh solution onto the coarse mesh.
     Hermes::Mixins::Loggable::Static::info("Calculating error estimate and exact error.");
-    OGProjection<double> ogProjection; ogProjection.project_global(&space, &ref_sln, &sln);
+    OGProjection<double> ogProjection; ogProjection.project_global(space, ref_sln, sln);
 
     // Calculate element errors and total error estimate.
-    Adapt<double> adaptivity(&space);
-    double err_est_rel = adaptivity.calc_err_est(&sln, &ref_sln) * 100;
-
     // Calculate exact error.
-    double err_exact_rel = Global<double>::calc_rel_error(&sln, &exact_sln, HERMES_H1_NORM) * 100;
+    errorCalculator.calculate_errors(sln, exact_sln, false);
+    double err_exact_rel = errorCalculator.get_total_error_squared() * 100;
+    errorCalculator.calculate_errors(sln, ref_sln);
+    double err_est_rel = errorCalculator.get_total_error_squared() * 100;
 
+    Adapt<double> adaptivity(space, &errorCalculator, &stoppingCriterion);
+
     cpu_time.tick();
     Hermes::Mixins::Loggable::Static::info("Error calculation: %g s", cpu_time.last());
 
     // Report results.
-    Hermes::Mixins::Loggable::Static::info("ndof_coarse: %d, ndof_fine: %d", space.get_num_dofs(), ref_space->get_num_dofs());
+    Hermes::Mixins::Loggable::Static::info("ndof_coarse: %d, ndof_fine: %d", space->get_num_dofs(), ref_space->get_num_dofs());
     Hermes::Mixins::Loggable::Static::info("err_est_rel: %g%%, err_exact_rel: %g%%", err_est_rel, err_exact_rel);
 
     // Time measurement.
     cpu_time.tick();
     double accum_time = cpu_time.accumulated();
 
     // View the coarse mesh solution and polynomial orders.
-    sview.show(&sln);
-    oview.show(&space);
+    sview.show(sln);
+    oview.show(space);
 
     // Add entry to DOF and CPU convergence graphs.
-    graph_dof_est.add_values(space.get_num_dofs(), err_est_rel);
+    graph_dof_est.add_values(space->get_num_dofs(), err_est_rel);
     graph_dof_est.save("conv_dof_est.dat");
     graph_cpu_est.add_values(accum_time, err_est_rel);
     graph_cpu_est.save("conv_cpu_est.dat");
-    graph_dof_exact.add_values(space.get_num_dofs(), err_exact_rel);
+    graph_dof_exact.add_values(space->get_num_dofs(), err_exact_rel);
     graph_dof_exact.save("conv_dof_exact.dat");
     graph_cpu_exact.add_values(accum_time, err_exact_rel);
     graph_cpu_exact.save("conv_cpu_exact.dat");
@@ -197,21 +194,17 @@ int main(int argc, char* argv[])
 
     // If err_est too large, adapt the mesh. The NDOF test must be here, so that the solution may be visualized
     // after ending due to this criterion.
-    if (err_exact_rel < ERR_STOP || space.get_num_dofs() >= NDOF_STOP) 
+    if (err_exact_rel < ERR_STOP || space->get_num_dofs() >= NDOF_STOP) 
       done = true;
     else
-      done = adaptivity.adapt(&selector, THRESHOLD, STRATEGY, MESH_REGULARITY);
+      done = adaptivity.adapt(&selector);
 
     cpu_time.tick();
     Hermes::Mixins::Loggable::Static::info("Adaptation: %g s", cpu_time.last());
 
     // Increase the counter of adaptivity steps.
     if (done == false)  
       as++;
-
-    if(done == false) 
-      delete ref_space->get_mesh();
-    delete ref_space;
   }
   while (done == false);
 

1d/moving-front/definitions.cpp

@@ -17,7 +17,7 @@ void CustomExactSolution::derivatives(double x, double y, double& dx, double& dy
   dy = 0;
 }
 
-Ord CustomExactSolution::ord(Ord x, Ord y) const 
+Ord CustomExactSolution::ord(double x, double y) const 
 {
   return Ord(20);
 }
@@ -43,23 +43,23 @@ CustomVectorFormVol::CustomVectorFormVol(int i,
   : VectorFormVol<double>(i), coeff(coeff), gt(gt)
 {
   this->set_area(area);
-  // If coeff is HERMES_ONE, initialize it to be constant 1.0.
-  if (coeff == HERMES_ONE) this->coeff = new Hermes::Hermes2DFunction<double>(1.0);
+  // If coeff is nullptr, initialize it to be constant 1.0.
+  if (coeff == nullptr) this->coeff = new Hermes::Hermes2DFunction<double>(1.0);
 }
 
 CustomVectorFormVol::CustomVectorFormVol(int i, 
 					 Hermes::Hermes2DFunction<double>* coeff, Hermes::vector<std::string> areas, GeomType gt)
   : VectorFormVol<double>(i), coeff(coeff), gt(gt)
 {
   this->set_areas(areas);
-  // If coeff is HERMES_ONE, initialize it to be constant 1.0.
-  if (coeff == HERMES_ONE) this->coeff = new Hermes::Hermes2DFunction<double>(1.0);
+  // If coeff is nullptr, initialize it to be constant 1.0.
+  if (coeff == nullptr) this->coeff = new Hermes::Hermes2DFunction<double>(1.0);
 }
 
 CustomVectorFormVol::~CustomVectorFormVol() 
 {
   // FIXME: Should be deleted here only if it was created here.
-  //if (coeff != HERMES_ONE) delete coeff;
+  //if (coeff != nullptr) delete coeff;
 };
 
 double CustomVectorFormVol::value(int n, double *wt, Func<double> *u_ext[], Func<double> *v,

1d/moving-front/definitions.h

@@ -1,5 +1,4 @@
 #include "hermes2d.h"
-#include "runge_kutta.h"
 
 using namespace Hermes;
 using namespace Hermes::Hermes2D;
@@ -9,15 +8,15 @@ using namespace Hermes::Hermes2D;
 class CustomExactSolution : public ExactSolutionScalar<double>
 {
 public:
-  CustomExactSolution(Mesh* mesh, double x0, double x1, double y0, double y1, 
+  CustomExactSolution(MeshSharedPtr mesh, double x0, double x1, double y0, double y1, 
                       double* t_ptr, double s, double c)
     : ExactSolutionScalar<double>(mesh), x0(x0), x1(x1), y0(y0), y1(y1), t_ptr(t_ptr), s(s), c(c) {};
 
   virtual double value(double x, double y) const;
 
   virtual void derivatives(double x, double y, double& dx, double& dy) const;
 
-  virtual Ord ord(Ord x, Ord y) const;
+  virtual Ord ord (double x, double y) const;
 
   double x0, x1, y0, y1, *t_ptr, s, c;
 };