Fix weak vs strong form (youpiee) (#30)

Commits:
* Documentation Fixes
* Add a check for normals ofedge and edgeInFront
* (Maybe) fix bad behaviour of the solution
N.B.: I still need to check the details

Author: ImperatorS79

Params/param.dat

@@ -1,6 +1,6 @@
-Gauss4
+Gauss3
 Lagrange
 RK1
 strong
-3
+300
 0.01

srcs/Mesh2D.cpp

@@ -97,6 +97,7 @@ static void loadElementProperties(std::map<int, ElementProperty>& meshElementPro
  * \param element The parent element.
  * \param nodesTagsEdge Node tags per edge of the element.
  * \param determinant1D Determinant associated with the edge of the element.
+ * \param nNodesElement Number of nodes of the parent element.
  */
 static void addEdge(Element2D& element, std::vector<int> nodesTagsEdge,
                     std::vector<double> determinant1D, unsigned int nNodesElement)
@@ -122,10 +123,9 @@ static void addEdge(Element2D& element, std::vector<int> nodesTagsEdge,
 
 
 /**
- * \brief Compute the outward normal of an edge
- * \param element the parent element
- * \param nodesTagsEdge Node tags per edge of the element
- * \param baryCenter  Barycenter of the parent element
+ * \brief Compute the outward normal of an edge.
+ * \param edge The edge of which the normal is computed.
+ * \param baryCenter  Barycenter of the parent element.
  */
 static void computeEdgeNormalCoord(Edge& edge,
                               const std::vector<double>& baryCenter)
@@ -170,10 +170,10 @@ static void computeEdgeNormalCoord(Edge& edge,
 }
 
 /**
- * \brief Compute the outward normal of an edge
+ * \brief Find the second position of the edge in the entity.
  * \param entity The parent entity.
  * \param currentEdge The edge of which you want to find the neighbor.
- * \param edgePos Index of the edge in his element.
+ * \param edgePos Index of the edge in its element.
  */
 static void findInFrontEdge(Entity2D& entity, Edge& currentEdge, unsigned int edgePos)
 {
@@ -213,6 +213,11 @@ static void findInFrontEdge(Entity2D& entity, Edge& currentEdge, unsigned int ed
     }
 }
 
+/**
+ * \brief Check if an edge lies on a boundary (and which one)
+ * \param nodesTagBoundaries Map which stores the tags of the nodes belonging to a certain boundary.
+ * \param edge The edge which we check if it is a boundary.
+ */
 static bool IsBounbdary(const std::map<std::string, std::vector<int>>& nodesTagBoundaries, Edge& edge)
 {
     for(std::pair<std::string, std::vector<int>> nodeTagBoundary : nodesTagBoundaries)
@@ -247,8 +252,9 @@ static bool IsBounbdary(const std::map<std::string, std::vector<int>>& nodesTagB
  * \param nGP1D Number of Gauss point for integration.
  * \param offsetInU Offset of the element in the u unknown vector.
  * \param nodesTagsPerEdge Node tags of the element, per edge.
- * \param intScheme Integration scheme for the basis functions evaluation.
- * \param basisFuncType The type of basis function you will use.
+ * \param elementBarycenter Element barycenter coordinate.
+ * \param nodesTagBoundaries Map which stores the tags of the nodes belonging to a certain boundary.
+ * \param element2DProperty Structure containing informations about a certain element type.
  */
 static void addElement(Entity2D& entity, int elementTag, int eleType2D,
                         int eleType1D, std::vector<double> jacobians2D,
@@ -257,9 +263,8 @@ static void addElement(Entity2D& entity, int elementTag, int eleType2D,
                         unsigned int nGP1D, unsigned int offsetInU,
                         const std::vector<int>& nodesTagsPerEdge,
                         const std::vector<double>& elementBarycenter,
-                        const std::string& intScheme,
-                        const std::string& basisFuncType,
-                        const std::map<std::string, std::vector<int>>& nodesTagBoundary)
+                        const std::map<std::string, std::vector<int>>& nodesTagBoundaries,
+                        const ElementProperty& element2DProperty)
 {
     Element2D element;
     element.elementTag = elementTag;
@@ -270,7 +275,6 @@ static void addElement(Entity2D& entity, int elementTag, int eleType2D,
 
     element.determinant2D = std::move(determinants2D);
     element.jacobian2D = std::move(jacobians2D);
-    unsigned int nNodesElement = nodesTagsPerEdge.size()/2; //Use elementProp map ?
 
     for(unsigned int i = 0 ; i < nodesTagsPerEdge.size()/2 ; ++i)
     {
@@ -279,16 +283,16 @@ static void addElement(Entity2D& entity, int elementTag, int eleType2D,
 
         std::vector<double> determinantsEdge1D(determinants1D.begin() + nGP1D*i, determinants1D.begin() + nGP1D*(i + 1));
 
-        addEdge(element, std::move(nodesTagsEdge), std::move(determinantsEdge1D), nNodesElement);
+        addEdge(element, std::move(nodesTagsEdge), std::move(determinantsEdge1D), element2DProperty.numNodes);
         if(entity.elements.size() != 0)
         {
-            if(!IsBounbdary(nodesTagBoundary, element.edges[i]))
+            if(!IsBounbdary(nodesTagBoundaries, element.edges[i]))
                 findInFrontEdge(entity, element.edges[i], i);
 
         }
         else
         {
-            IsBounbdary(nodesTagBoundary, element.edges[i]);
+            IsBounbdary(nodesTagBoundaries, element.edges[i]);
         }
 
         computeEdgeNormalCoord(element.edges[i], elementBarycenter);
@@ -396,8 +400,9 @@ static void addEntity(Mesh2D& mesh, const std::pair<int, int>& entityHandle, uns
                         std::move(determinantElement1D),
                         nGP1D, currentOffset,
                         nodesTagPerEdgeElement,
-                        elementBarycenter, intScheme, basisFuncType,
-                        mesh.nodesTagBoundary);
+                        elementBarycenter,
+                        mesh.nodesTagBoundary,
+                        mesh.elementProperties2D[eleType2D]);
 
             currentOffset += nodesTagPerEdgeElement.size()/2;
         }

srcs/Mesh2D.hpp

@@ -29,7 +29,7 @@ struct Edge
     std::string bcName; /**< Name of the physical group of the boundary condition
                              in which the edge is in (if any)*/
 
-    std::vector<unsigned int> offsetInElm;
+    std::vector<unsigned int> offsetInElm;  /**< Offset of the edge nodes in the element*/
 };
 
 /**
@@ -42,7 +42,7 @@ struct Element2D
     int elementType2D;                  /**< 2D type of the element*/
     int elementType1D;                  /**< 1D type of the element*/ //Store it only once
 
-    unsigned int offsetInU;
+    unsigned int offsetInU;             /**< Offset of the element in the unknowns vector*/
 
     std::vector<double> determinant2D;  /**< Determinant of the variable change for Gauss integration,
                                              evaluated at each Gauss point*/
@@ -103,8 +103,7 @@ struct Mesh2D
     std::map<int, ElementProperty> elementProperties2D; /**< Store element properties for each 2D element type */
     std::map<int, ElementProperty> elementProperties1D; /**< Store element properties for each 1D element type */
 
-    std::map<std::string, std::vector<int>> nodesTagBoundary;
-    //std::map<std::string, std::vector<double>> coordNodesBoundary;
+    std::map<std::string, std::vector<int>> nodesTagBoundary;/**< Tags of the nodes per BC */
 
     std::vector<Entity2D> entities; /**< List of entities inside the mesh */
 };
@@ -113,13 +112,15 @@ struct Mesh2D
 /**
  * \brief Get the number of nodes (i.e. of unknowns) given a mesh.
  * \param mesh2D The structure that contains the mesh.
+ * \return The number of unknowns inside a mesh.
  */
 unsigned int getNumNodes(const Mesh2D& mesh2D);
 
 
 /**
  * \brief Get the tags of nodes (i.e. of unknowns) given a mesh.
  * \param mesh2D The structure that contains the mesh.
+ * \return Vector containing all the tags of the nodes inside a mesh.
  */
 std::vector<int> getTags(const Mesh2D& mesh2D);
 

srcs/buildFlux.cpp

@@ -7,25 +7,23 @@ void flux(Eigen::VectorXd& fx, Eigen::VectorXd& fy, double& C,
 			const Eigen::VectorXd& u)
 {
 	// first basic flux: simple transport
-	std::vector<double> a;
-	a.push_back(1.0);
-	a.push_back(0.0);
+	double ax = 1.0;
+	double ay = 0.0;
 
-	fx = a[0]*u;
-	fy = a[1]*u;
+	fx = ax*u;
+	fy = ay*u;
 
-	C = sqrt(a[0]*a[0] + a[1]*a[1]);
+	C = sqrt(ax*ax + ay*ay);
 }
 
 void flux(double& fx, double& fy, double u)
 {
 	// first basic flux: simple transport
-	std::vector<double> a;
-	a.push_back(1.0);
-	a.push_back(0.0);
+	double ax = 1.0;
+	double ay = 0.0;
 
-	fx = a[0]*u;
-	fy = a[1]*u;
+	fx = ax*u;
+	fy = ay*u;
 }
 
 
@@ -34,7 +32,7 @@ double valueAtBC(const std::string& bcName,
 {
 	if(!bcName.compare("BC_Left"))
 	{
-		return 0.01;//sin(20*t);
+		return exp(-(t-2)*(t-2)/0.1);
 	}
 	else
 	{
@@ -45,19 +43,10 @@ double valueAtBC(const std::string& bcName,
 
 bool buildFlux(const Mesh2D& mesh, Eigen::VectorXd& I, const Eigen::VectorXd& u,
 	const Eigen::VectorXd& fx, const Eigen::VectorXd& fy, const double& C,
-	const std::string& typeForm, unsigned int numNodes, double t)
+	double factor, unsigned int numNodes, double t)
 {
 
 	// the type of form is stored in factor
-	float factor;
-	if(!typeForm.compare("strong")) 	factor = -1.0;
-	else if(!typeForm.compare("weak"))	factor = +1.0;
-	else
-	{
-		std::cerr 	<< "The form  " << typeForm  << "does not exist !"
-					<< std::endl;
-		return false;
-	}
 
 	// loop over the entities
 	for(unsigned int ent = 0 ; ent < mesh.entities.size() ; ent++)
@@ -136,11 +125,11 @@ bool buildFlux(const Mesh2D& mesh, Eigen::VectorXd& I, const Eigen::VectorXd& u,
 							edge.nodeCoordinate[j].second,
 							0.0, t);
 
-						flux(fxAtBC, fyAtBC, uAtBC);
+                        flux(fxAtBC, fyAtBC, uAtBC);
 
-						gx[j] = -(factor*fx[indexJ] + fxAtBC)/2
+						gx[edge.offsetInElm[j]] += -(factor*fx[indexJ] + fxAtBC)/2
 							- C*element.edges[s].normal.first*(u[indexJ] - uAtBC)/2;
-						gy[j] = -(factor*fy[indexJ] + fyAtBC)/2
+						gy[edge.offsetInElm[j]] += -(factor*fy[indexJ] + fyAtBC)/2
 							- C*element.edges[s].normal.second*(u[indexJ] - uAtBC)/2;
 					}
 					else
@@ -154,9 +143,9 @@ bool buildFlux(const Mesh2D& mesh, Eigen::VectorXd& I, const Eigen::VectorXd& u,
                        					.edges[edge.edgeInFront.second]
                        					.offsetInElm[edge.nodeIndexEdgeInFront[j]];
 
-						gx[j] = -(factor*fx[indexJ] + fx[indexFrontJ])/2
+						gx[edge.offsetInElm[j]] += -(factor*fx[indexJ] + fx[indexFrontJ])/2
 							- C*element.edges[s].normal.first*(u[indexJ] - u[indexFrontJ])/2;
-						gy[j] = -(factor*fy[indexJ] + fy[indexFrontJ])/2
+						gy[edge.offsetInElm[j]] += -(factor*fy[indexJ] + fy[indexFrontJ])/2
 							- C*element.edges[s].normal.second*(u[indexJ] - u[indexFrontJ])/2;
 						std::cout << "==========>" << gx[j] << std::endl;
 					}

srcs/buildFlux.hpp

@@ -12,6 +12,6 @@ void flux(double& fx, double& fy, double u);
 
 bool buildFlux(const Mesh2D& mesh, Eigen::VectorXd& I, const Eigen::VectorXd& u,
 	const Eigen::VectorXd& fx, const Eigen::VectorXd& fy, const double& C,
-	const std::string& typeForm, unsigned int numNodes, double t);
+	double factor, unsigned int numNodes, double t);
 
 #endif /* buildFlux_hpp */

srcs/displayMesh.cpp

@@ -80,6 +80,12 @@ void displayMesh(const Mesh2D& mesh)
                                 << edges[k].edgeInFront.first<<", "
                                 << "edge " << edges[k].edgeInFront.second<<", "
                                 << std::endl; //[TO DO]: put again inverted
+                    if((edges[k].normal.first != -entitity.elements[edges[k].edgeInFront.first].edges[edges[k].edgeInFront.second].normal.first)
+                       || (edges[k].normal.second != -entitity.elements[edges[k].edgeInFront.first].edges[edges[k].edgeInFront.second].normal.second))
+                    {
+                        std::cerr<<"Bug in the normal of that edge !"<<std::endl;
+                        return;
+                    }
                 }
                 else
                 {