Add a check for Lagrange polynomials order

And replace x, y, z by a pos vector for ibvFunction

Author: ImperatorS79

srcs/mesh/Mesh.cpp

@@ -389,6 +389,7 @@ static void addElement(Entity& entity, int elementTag, int eleTypeHD,
  * \param currentOffset Offset in u unknown vector.
  * \param intScheme Integration scheme for the basis functions evaluation.
  * \param basisFuncType The type of basis function you will use.
+ * \return true if then entity was added flawlessly, false otherwise.
  */
 static bool addEntity(Mesh& mesh, int entityTag, unsigned int& currentOffset,
                       const std::string& intScheme, const std::string& basisFuncType)
@@ -424,6 +425,14 @@ static bool addEntity(Mesh& mesh, int entityTag, unsigned int& currentOffset,
         unsigned int numNodes = mesh.elementProperties[eleTypeHD].numNodes;
         unsigned int order = mesh.elementProperties[eleTypeHD].order;
 
+        if(intScheme == "Lagrange" && order > 7)
+        {
+            std::cerr << "Lagrange polynomials are not stable with an order superior to 7"
+                      << std::endl;
+
+            return false;
+        }
+
         //Get the elements barycenters (for normal computation)
         std::vector<double> baryCenters;
         gmsh::model::mesh::getBarycenters(eleTypeHD, entityTag, false, true, baryCenters);

srcs/params/Params.cpp

@@ -115,8 +115,9 @@ static bool handleBoundaryCondition(std::ifstream& paramFile, SolverParams& solv
         return false;
     }
 
-    std::cout << "Inital Condition present and " << nBC
-              << " boundary conditions present in file " << fileName << std::endl;
+    std::cout << "Initial condition present and " << nBC
+              << " boundary conditions present in "<< std::endl
+              << "file " << fileName << std::endl;
 
     return true;
 }

srcs/params/ibvFunction.cpp

@@ -2,43 +2,44 @@
 #include <cassert>
 #include "ibvFunction.hpp"
 
-double sinus(double x, double y, double z, double u,
+double sinus(const std::vector<double>& pos, double u,
              double t, const std::vector<double>& coeffs)
 {
     assert(coeffs.size() == 3);
 
     return coeffs[0]*sin(2*M_PI*coeffs[1]*t + coeffs[2]);
 }
 
-double gaussian(double x, double y, double z, double u,
+double gaussian(const std::vector<double>& pos, double u,
                 double t, const std::vector<double>& coeffs)
 {
     assert(coeffs.size() == 3);
 
     return coeffs[0]*exp(-(t-coeffs[1])*(t-coeffs[1])/(2*coeffs[2]));
 }
 
-double constant(double x, double y, double z, double u,
+double constant(const std::vector<double>& pos, double u,
                 double t, const std::vector<double>& coeffs)
 {
     assert(coeffs.size() == 1);
 
     return coeffs[0];
 }
 
-double constantNeumann(double x, double y, double z, double u,
+double constantNeumann(const std::vector<double>& pos, double u,
                 double t, const std::vector<double>& coeffs)
 {
     return u;
 }
 
-double gaussian2D(double x, double y, double z, double u,
+double gaussian2D(const std::vector<double>& pos, double u,
                   double t, const std::vector<double>& coeffs)
 {
     assert(coeffs.size() == 5);
+    assert(pos.size() == 3);
 
-    double X = (x-coeffs[1])*(x-coeffs[1])/(2*coeffs[2]);
-    double Y = (y-coeffs[3])*(y-coeffs[3])/(2*coeffs[4]);
+    double X = (pos[0]-coeffs[1])*(pos[0]-coeffs[1])/(2*coeffs[2]);
+    double Y = (pos[1]-coeffs[3])*(pos[1]-coeffs[3])/(2*coeffs[4]);
 
     return coeffs[0]*exp(-(X+Y));
 }

srcs/params/ibvFunction.hpp

@@ -12,78 +12,68 @@
 struct ibc
 {
     std::vector<double> coefficients;   /**< Coefficient for the mathematical function */
-    std::function<double(double x, double y, double z,
+    std::function<double(const std::vector<double>& pos,
                          double u, double t,
                          std::vector<double> coeffs)> ibcFunc; /**< Pointer to the mathematical function */
 };
 
 
 /**
  * \brief Compute a simple A*sin(2*pi*nu*t + phi)
- * \param x x coordinate.
- * \param y y coordinate.
- * \param z z coordinate.
+ * \param pos Node position.
  * \param u The current solution.
  * \param t Current time.
  * \param coeffs Coefficient for the sinus. coeffs[0] = A, coeffs[1] = nu
  * coeffs[2] = phi.
  * \return Value of the function at (x, y , z, t).
  */
-double sinus(double x, double y, double z, double u,
+double sinus(const std::vector<double>& pos, double u,
              double t, const std::vector<double>& coeffs);
 
 /**
  * \brief Compute a simple A*exp(-(t-t_peak)^2/(2*var))
- * \param x x coordinate.
- * \param y y coordinate.
- * \param z z coordinate.
+ * \param pos Node position.
  * \param u The current solution.
  * \param t Current time.
  * \param coeffs Coefficient for the gaussian. coeffs[0] = A, coeffs[1] = t_peak
  * coeffs[2] = var.
  * \return Value of the function at (x, y , z, t).
  */
-double gaussian(double x, double y, double z, double u,
+double gaussian(const std::vector<double>& pos, double u,
                 double t, const std::vector<double>& coeffs);
 
 /**
  * \brief Compute a constant
- * \param x x coordinate.
- * \param y y coordinate.
- * \param z z coordinate.
+ * \param pos Node position.
  * \param u The current solution.
  * \param t Current time.
  * \param coeffs Coefficient for the constant. coeffs[0] = constant.
  * \return Value of the function at (x, y , z, t).
  */
-double constant(double x, double y, double z, double u,
+double constant(const std::vector<double>& pos, double u,
                 double t, const std::vector<double>& coeffs);
 
 /**
  * \brief Compute a Von Neumann constant value
- * \param x x coordinate.
- * \param y y coordinate.
- * \param z z coordinate.
+ * \param pos Node position.
  * \param u The current solution.
  * \param t Current time.
  * \param coeffs Coefficient (not used here).
  * \return u.
  */
-double constantNeumann(double x, double y, double z, double u,
+double constantNeumann(const std::vector<double>& pos, double u,
                        double t, const std::vector<double>& coeffs);
 
 /**
  * \brief Compute a simple A*exp(-(x-x0)^2/(2*var_y) -(y-y0)^2/(2*var_y))
- * \param x x coordinate.
- * \param y y coordinate.
- * \param z z coordinate.
+ * \param pos Node position.
  * \param u The current solution.
  * \param t Current time.
  * \param coeffs Coefficient for the gaussian. coeffs[0] = A, coeffs[1] = x_0
  * coeffs[2] = var_x, coeffs[3]=y_0, coeffs[4]=var_y.
  * \return Value of the function at (x, y , z, t).
  */
-double gaussian2D(double x, double y, double z, double u,
+double gaussian2D(const std::vector<double>& pos, double u,
                   double t, const std::vector<double>& coeffs);
 
 #endif // bcFunction_hpp_included

srcs/solver/buildFlux.cpp

@@ -83,9 +83,8 @@ void buildFlux(const Mesh& mesh, Eigen::VectorXd& I, const Eigen::VectorXd& u,
 						ibc boundary = boundaries.at(edge.bcName);
 
 						// node "in front" (at boundary condition)
-						uAtBC = boundary.ibcFunc(edge.nodeCoordinate[j][0],
-							edge.nodeCoordinate[j][1],
-							0.0, u[indexJ], t, boundary.coefficients);
+						uAtBC = boundary.ibcFunc(edge.nodeCoordinate[j],
+                               u[indexJ], t, boundary.coefficients);
 
 						// physical flux "in front" (at boundary condition)
                         flux(fxAtBC, fyAtBC, uAtBC, fluxCoeffs);

srcs/solver/timeInteg.cpp

@@ -145,10 +145,8 @@ bool timeInteg(const Mesh& mesh, const SolverParams& solverParams,
         {
             for(unsigned int n = 0 ; n < element.nodeTags.size() ; ++n)
             {
-                double x = element.nodesCoord[n][0];
-                double y = element.nodesCoord[n][1];
                 u(element.offsetInU + n) =
-                solverParams.initCondition.ibcFunc(x, y, 0, 0, 0,
+                solverParams.initCondition.ibcFunc(element.nodesCoord[n], 0, 0,
                                                    solverParams.initCondition.coefficients);
             }
         }