diff --git a/srcs/BEM/boundaryElements.cpp b/srcs/BEM/boundaryElements.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bdfa2d7408deecb5ec785e2be9e745c6d5fc9eb2
--- /dev/null
+++ b/srcs/BEM/boundaryElements.cpp
@@ -0,0 +1,311 @@
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+#include <sstream>
+#include <Eigen/Dense>
+#include <cmath>
+#include <iomanip>
+#include <stdio.h>
+#include <chrono>
+#include "functionsBEM.hpp"
+#include <omp.h>
+#include <algorithm>
+#include <list>
+
+// Fills the "boundaryConditions" map associating to each physical group (a string) the type of boundary condition ('true':
+// dirichlet, 'false': neumann) and its value (if any). "domainName" is the name of the current BEM domain.
+void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
+ const std::string domainName)
+{
+ // Extracts all the boundary conditions from "ONELAB".
+ std::vector<std::string> bcKeys;
+ gmsh::onelab::getNames(bcKeys, "Boundary Conditions.+");
+
+ for(auto &key : bcKeys)
+ {
+ std::vector<double> values;
+ gmsh::onelab::getNumber(key, values);
+
+ std::stringstream ss(key);
+ std::vector<std::string> words;
+ std::string word;
+ while(getline(ss, word, '/'))
+ words.push_back(word);
+
+ // Only focus on the right BEM domain.
+ if(words.size() == 4 && words[2].compare(domainName) == 0)
+ {
+ if(words[3].compare("dirichlet") == 0)
+ boundaryConditions[words[1]] = {true, values[0]};
+ else if(words[3].compare("neumann") == 0)
+ boundaryConditions[words[1]] = {false, values[0]};
+ }
+ }
+}
+
+// Fills the "elementVector" vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the
+// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also fills the
+// "physicalGroups" map associating to each (physical group) name, its tag and dimension. "domainName" is the name of the current
+// BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement to a serie of node tags. Returns the number of
+// 'dirichlet elements' (useful for the construction/resolution of the system).
+std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
+ std::map<std::string, std::pair<int, int>> &physicalGroups,
+ std::map<int, std::pair<double, double>> &nodalDisplacements,
+ const std::string domainName)
+{
+ std::size_t nbDirichletElements = 0; // Will be incremented later
+
+ std::map<std::string, std::pair<bool, double>> boundaryConditions;
+ getBC(boundaryConditions, domainName);
+
+ // Gets the dim and the tags of the physical groups.
+ gmsh::vectorpair physicalGroupsDimTags;
+ gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
+
+ for(std::size_t i = 0; i < physicalGroupsDimTags.size(); i++)
+ {
+ // Fills the "physicalGroups" map that associates the name to the {dim, tag}.
+ int physicalGroupDim = physicalGroupsDimTags[i].first;
+ int physicalGroupTag = physicalGroupsDimTags[i].second;
+ std::string physicalGroupName;
+ gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
+ physicalGroups[physicalGroupName] = {physicalGroupDim, physicalGroupTag};
+
+ gmsh::model::getPhysicalName(physicalGroupDim, physicalGroupTag, physicalGroupName);
+
+ // If no boundary condition on a given physical curve, don't do anything.
+ if(boundaryConditions.find(physicalGroupName) == boundaryConditions.end())
+ continue;
+
+ // Retrieves the type and the value of the boundary condition of this physical group.
+ bool bcType = boundaryConditions[physicalGroupName].first;
+ double bcValue = boundaryConditions[physicalGroupName].second;
+
+ // Loops over the tags of the entities and the types of the elements for retrieving the relevant information (element
+ // tags and nodes tags) in gmsh.
+ std::vector<int> entityTags;
+ gmsh::model::getEntitiesForPhysicalGroup(physicalGroupDim, physicalGroupTag, entityTags);
+ for (std::size_t i = 0; i < entityTags.size(); i++)
+ {
+ std::vector<int> elementTypes = {1, 8}; // Hardcoded types of 2/3-nodes 1D elements.
+ for(std::size_t j = 0; j < elementTypes.size(); j++)
+ {
+ int nbNodesPerElement = 2;
+ if(elementTypes[j] == 8)
+ nbNodesPerElement = 3;
+
+ // Get the element tags and the node tags of each element.
+ std::vector<std::size_t> elementTags;
+ std::vector<std::size_t> elementNodeTags;
+ gmsh::model::mesh::getElementsByType(elementTypes[j], elementTags, elementNodeTags, entityTags[i]);
+
+ // Get the coordinates of each node (extracted above).
+ std::vector<std::size_t> nodeTags;
+ std::vector<double> nodeCoords;
+ std::vector<double> parametricNodeCoords;
+ gmsh::model::mesh::getNodesByElementType(elementTypes[j], nodeTags, nodeCoords, parametricNodeCoords, entityTags[i], false);
+
+ // Constructs a map associating the index of each node tag in the "nodeCoords" vector.
+ std::map<std::size_t, std::size_t> nodeIndices;
+ for(std::size_t k = 0; k < nodeTags.size(); k++)
+ nodeIndices[nodeTags[k]] = k;
+
+ // Computes each element data and fills the "elementVector" vector.
+ #pragma omp parallel shared(nbDirichletElements, elementTags, nodeTags, nodalDisplacements, nodeIndices, nodeCoords, bcType, bcValue)
+ {
+ std::vector<elementStruct> localDirichletElementVector;
+ std::vector<elementStruct> localNeumannElementVector;
+
+ // Computes the element data
+ #pragma omp for schedule(guided)
+ for(std::size_t el = 0; el < elementTags.size(); el++)
+ {
+ std::size_t elTag = elementTags[el];
+
+ elementStruct element;
+ element.tag = elTag;
+ element.nodeTags[0] = nodeTags[nbNodesPerElement*el];
+ element.nodeTags[1] = nodeTags[nbNodesPerElement*el + 1];
+
+ if(nbNodesPerElement == 3)
+ element.midNodeTag = nodeTags[nbNodesPerElement*el + 2];
+ else
+ element.midNodeTag = 0; // "midNodeTag" set to '0' (default value) if no mid node on the element.
+
+ std::pair<double, double> nodalDisp1 = {0.0, 0.0}, nodalDisp2 = {0.0, 0.0};
+ if(nodalDisplacements.count(element.nodeTags[0]))
+ nodalDisp1 = nodalDisplacements[element.nodeTags[0]];
+ if(nodalDisplacements.count(element.nodeTags[1]))
+ nodalDisp2 = nodalDisplacements[element.nodeTags[1]];
+ computeGeometricParam(element, nodeIndices, nodeCoords, nodalDisp1, nodalDisp2);
+
+ // Pushes the element at the back of lists (private for each thread)/
+ if(bcType)
+ {
+ element.bcValue[0] = bcValue;
+ localDirichletElementVector.push_back(element);
+ }
+ else
+ {
+ element.bcValue[1] = bcValue;
+ localNeumannElementVector.push_back(element);
+ }
+ }
+
+ // Add the elements to the "elementVector" vector. It is important that the 'dirichlet elements' are first in
+ // the vector.
+ #pragma omp critical
+ {
+ elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(), localDirichletElementVector.end());
+ elementVector.insert(elementVector.end(), localNeumannElementVector.begin(), localNeumannElementVector.end());
+ nbDirichletElements += localDirichletElementVector.size();
+ }
+ }
+ }
+ }
+ }
+
+ return nbDirichletElements;
+}
+
+// Computes the useful data of a given element and stores it in "element". "nodeIndinces" is a map associating the index of any
+// node tag for accessing its coordinates in the "nodeCoords" vector. "nodalDisp1" and "nodalDisp2" are pairs representing the
+// potential displacement of the first and second node (respectivelly) of the element. If an element has 3 nodes, the first and
+// second nodes are considered to be the extremities of that element.
+void computeGeometricParam(elementStruct &element,
+ std::map<std::size_t, std::size_t> &nodeIndices,
+ const std::vector<double> &nodeCoords,
+ const std::pair<double, double> &nodalDisp1,
+ const std::pair<double, double> &nodalDisp2)
+{
+ std::size_t i1Tag = element.nodeTags[0];
+ std::size_t i2Tag = element.nodeTags[1];
+
+ std::size_t i1Index = nodeIndices[i1Tag];
+ std::size_t i2Index = nodeIndices[i2Tag];
+
+ // Get the coordinates of the nodes from their tag.
+ element.x1 = nodeCoords[3*i1Index] + nodalDisp1.first;
+ element.x2 = nodeCoords[3*i2Index] + nodalDisp2.first;
+ element.y1 = nodeCoords[3*i1Index + 1] + nodalDisp1.second;
+ element.y2 = nodeCoords[3*i2Index + 1] + nodalDisp2.second;
+
+ element.midX = (element.x1 + element.x2) / 2;
+ element.midY = (element.y1 + element.y2) / 2;
+ element.deltaX = element.x2 - element.x1;
+ element.deltaY = element.y2 - element.y1;
+ element.length = sqrt(element.deltaX * element.deltaX + element.deltaY * element.deltaY);
+}
+
+// Fills the "A", "B", "c" matrices (/vector) that represent the linear system. "elementVector" is a vector storing the boundary
+// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "nbDirichletElements" is the number of
+// 'dirichlet elements'. "solutionType" indicates whether the elements of the matrices are computed analytically ('true') or
+// numerically ('false'). "integrationType" is the type of numerical integration used (in case of numerical computation).
+void fillSystem(Eigen::MatrixXd &A,
+ Eigen::MatrixXd &B,
+ Eigen::MatrixXd &c,
+ const std::vector<elementStruct> &elementVector,
+ const std::size_t nbDirichletElements,
+ const std::string integrationType,
+ const bool solutionType)
+{
+ double h;
+ double g;
+
+ // Double loop on the elements for computing each elements of the linear system.
+ #pragma omp parallel for collapse(2) schedule(guided) shared(elementVector, A, B, c) private(h, g)
+ for(std::size_t i = 0; i < elementVector.size(); i++)
+ {
+ for(std::size_t j = 0; j < elementVector.size(); j++)
+ {
+ // Will be computed and stored in A and B.
+ h = 0;
+ g = 0;
+
+ if(i == j) // Formula for diagonal elements.
+ {
+ h = -0.5;
+ g = (elementVector[i].length/(2*M_PI)) * (log(elementVector[i].length/2)-1);
+ }
+ else // General case
+ {
+ h = computeAnalyticalH(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
+ if(!solutionType)
+ g = computeNumericalG(elementVector[i].midX, elementVector[i].midY, elementVector[j], integrationType);
+ else
+ g = computeAnalyticalG(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
+ }
+
+ // Storage of h and g at the right place (in matrices A and B and vector c).
+ if(j < nbDirichletElements)
+ {
+ A(i,j) = g;
+ B(i,j) = h;
+ c(j) = elementVector[j].bcValue[0];
+ }
+ else
+ {
+ A(i,j) = -h;
+ B(i,j) = -g;
+ c(j) = elementVector[j].bcValue[1];
+ }
+ }
+ }
+}
+
+// Fills the "electrostaticPressure" map associating each (boundary) element tag with its electrostatic pressure. Also sets the
+// unknown of each element (either the potential of the electric field) of the "elementVector" vector storing the boundary
+// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "x" is the solution of the linear
+// system. "nbDirichletElements" is the number of 'dirichlet elements', "domainName" is the name of the considered BEM domain.
+void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
+ std::vector<elementStruct> &elementVector,
+ const Eigen::MatrixXd &x,
+ const std::size_t nbDirichletElements,
+ const std::string domainName)
+{
+ // Retrieves the dielectric permittivity of the BEM domain. If not set in the .geo file, set to the dielectric permittivity
+ // of vacuum (with a warning).
+ std::vector<std::string> keys;
+ gmsh::onelab::getNames(keys, "(Materials).+");
+
+ double epsilon = 8.8541878128e-12; // dielectric permittivity
+ bool setEpsilon = false;
+ for (auto &key : keys)
+ {
+ std::vector<double> value;
+ gmsh::onelab::getNumber(key, value);
+
+ // Expected key structure is "type/group_name/field"
+ std::stringstream ss(key);
+ std::vector<std::string> words;
+ std::string word;
+ while(std::getline(ss, word, '/')) // read string until '/'
+ words.push_back(word);
+ if(words.size()==3)
+ {
+ if(words[2].compare("Epsilon") == 0 && words[1].compare(domainName) == 0)
+ {
+ epsilon = value[0];
+ setEpsilon = true;
+ }
+ }
+ }
+ if(!setEpsilon)
+ std::cout << "\tWarning: You forgot to set the dielectric permittivity of the material for the domain '" << domainName << "' in the .geo file and it thus has been fixed to the dielectric permittivity of the vacuum: " << epsilon << " [F/m].\n";
+
+ #pragma omp parallel for schedule(guided)
+ for(size_t i = 0; i < elementVector.size(); i++)
+ {
+ // Stores the unknown value (potential or electric field) in the appropriate elements.
+ if(i < nbDirichletElements)
+ elementVector[i].bcValue[1] = x(i);
+ else
+ elementVector[i].bcValue[0] = x(i);
+
+ // Computes the electrostatic pressure and fills the map.
+ double pressure = epsilon * elementVector[i].bcValue[1] * elementVector[i].bcValue[1] / 2; // don't forget factor 2
+
+ #pragma omp critical
+ electrostaticPressure[elementVector[i].tag] = pressure;
+ }
+}
\ No newline at end of file
diff --git a/srcs/BEM/computations.cpp b/srcs/BEM/computations.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c35b66d0c76ecebd80aff3deae8237ecc8571920
--- /dev/null
+++ b/srcs/BEM/computations.cpp
@@ -0,0 +1,330 @@
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+#include <sstream>
+#include <Eigen/Dense>
+#include <cmath>
+#include <iomanip>
+#include <stdio.h>
+#include <chrono>
+#include "functionsBEM.hpp"
+#include <omp.h>
+#include <algorithm>
+#include <list>
+
+double computeAnalyticalH(const double x,
+ const double y,
+ const elementStruct &element)
+{
+ double deltaX1 = element.x1 - x;
+ double deltaX2 = element.x2 - x;
+ double deltaY1 = element.y1 - y;
+ double deltaY2 = element.y2 - y;
+ double cosAlpha = (deltaX1*deltaX2 + deltaY1*deltaY2) / (sqrt((pow(deltaX1,2.0) + pow(deltaY1,2.0)) * (pow(deltaX2,2.0) + pow(deltaY2,2.0))));
+
+ if(cosAlpha < -1.0)
+ {
+ cosAlpha = -1.0;
+ }
+ else if(cosAlpha > 1.0)
+ {
+ cosAlpha = 1.0;
+ }
+
+ double alpha = acos(cosAlpha);
+ double crossProduct = deltaX1*deltaY2 - deltaY1*deltaX2;
+
+ if(crossProduct < 0)
+ {
+ alpha = -alpha;
+ }
+
+ return alpha / (2*M_PI);
+}
+double computeAnalyticalG(const double x,
+ const double y,
+ const elementStruct &element)
+{
+ double lj = element.length;
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+ double Pi = M_PI;
+
+ double g_analytique = (lj*(-2*(pow(xj - xjp1,2) + pow(yj - yjp1,2)) + 2*(-(xj*yi) + xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*yjp1)* atan(((xi - xjp1)*(xj - xjp1) + (yi - yjp1)*(yj - yjp1))/(-(xj*yi) + xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*yjp1)) + 2*(-(xjp1*yi) - xi*yj + xjp1*yj + xj*(yi - yjp1) + xi*yjp1)* atan(((xi - xj)*(xj - xjp1) + (yi - yj)*(yj - yjp1))/ (xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1))) + (-((xi - xj)*(xj - xjp1)) - (yi - yj)*(yj - yjp1))* log(pow(xi - xj,2) + pow(yi - yj,2)) + ((xi - xjp1)*(xj - xjp1) + (yi - yjp1)*(yj - yjp1))* log(pow(xi - xjp1,2) + pow(yi - yjp1,2))))/(4.*Pi*(pow(xj - xjp1,2) + pow(yj - yjp1,2)));
+
+ return g_analytique;
+}
+double computeAnalyticalGradHx(const double x,
+ const double y,
+ const elementStruct &element)
+{
+ double lj = element.length;
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+ double Pi = M_PI;
+
+ double dx = xjp1-xj;
+ double dy = yjp1-yj;
+ double nx = dy;
+ double ny = -dx;
+ double norm = sqrt(nx*nx+ny*ny);
+ nx = nx / norm;
+ ny = ny / norm;
+
+ double gradHx_analytique = (lj*(nx*(pow(xi,2) + xj*xjp1 - xi*(xj + xjp1) - (yi - yj)*(yi - yjp1)) -
+ ny*(xj*yi + xjp1*yi - xjp1*yj - xj*yjp1 + xi*(-2*yi + yj + yjp1))))/
+ (2.*Pi*(pow(xi - xj,2) + pow(yi - yj,2))*
+ (pow(xi - xjp1,2) + pow(yi - yjp1,2)));
+
+ return gradHx_analytique;
+}
+double computeAnalyticalGradHy(const double x,
+ const double y,
+ const elementStruct &element)
+{
+ double lj = element.length;
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+ double Pi = M_PI;
+
+ double dx = xjp1-xj;
+ double dy = yjp1-yj;
+ double nx = dy;
+ double ny = -dx;
+ double norm = sqrt(nx*nx+ny*ny);
+ nx = nx / norm;
+ ny = ny / norm;
+
+ double gradHy_analytique = -(lj*(ny*(pow(xi,2) + xj*xjp1 - xi*(xj + xjp1) - (yi - yj)*(yi - yjp1)) +
+ nx*(xj*yi + xjp1*yi - xjp1*yj - xj*yjp1 + xi*(-2*yi + yj + yjp1))))/
+ (2.*Pi*(pow(xi - xj,2) + pow(yi - yj,2))*
+ (pow(xi - xjp1,2) + pow(yi - yjp1,2)));
+
+ return gradHy_analytique;
+}
+double computeAnalyticalGradGx(const double x,
+ const double y,
+ const elementStruct &element)
+{
+ double lj = element.length;
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+ double Pi = M_PI;
+
+ double gradGx_analytique = (lj*(-2*(yj - yjp1)*(atan(((xi - xj)*(xj - xjp1) + (yi - yj)*(yj - yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1))) +
+ atan(((xi - xjp1)*(-xj + xjp1) + (yi - yjp1)*(-yj + yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1)))) -
+ (xj - xjp1)*(log(pow(xi - xj,2) + pow(yi - yj,2)) - log(pow(xi - xjp1,2) + pow(yi - yjp1,2)))))/(4.*Pi*(pow(xj - xjp1,2) + pow(yj - yjp1,2)));
+
+ return gradGx_analytique;
+}
+double computeAnalyticalGradGy(const double x,
+ const double y,
+ const elementStruct &element)
+{
+ double lj = element.length;
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+ double Pi = M_PI;
+
+ double gradGy_analytique = (lj*(2*(xj - xjp1)*(atan(((xi - xj)*(xj - xjp1) + (yi - yj)*(yj - yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1))) +
+ atan(((xi - xjp1)*(-xj + xjp1) + (yi - yjp1)*(-yj + yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1)))) -
+ (yj - yjp1)*(log(pow(xi - xj,2) + pow(yi - yj,2)) - log(pow(xi - xjp1,2) + pow(yi - yjp1,2)))))/(4.*Pi*(pow(xj - xjp1,2) + pow(yj - yjp1,2)));
+ return gradGy_analytique;
+}
+
+
+double computeNumericalG(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType)
+{
+ // Computation of the integration points (parametric) and weights.
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
+
+ // Computation of the sum (of integrand evaluated at integration point * weight)
+ double gSum = 0;
+ for(size_t k = 0; k < weights.size(); k++)
+ {
+ // Parametric coordinate
+ double xi = localCoord[3*k];
+ // Global coordinates
+ double gX = element.midX + element.deltaX*xi/2;
+ double gY = element.midY + element.deltaY*xi/2;
+
+ // Distance between mid point of element 'i' and integration point.
+ double r = sqrt(pow(x - gX, 2.0) + pow(y - gY, 2.0));
+ gSum += log(r) * weights[k];
+ }
+
+ double g = gSum * element.length / (4*M_PI);
+
+ return g;
+}
+double computeNumericalGradHx(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType)
+{
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+
+ double dx = xjp1-xj;
+ double dy = yjp1-yj;
+ double nx = dy;
+ double ny = -dx;
+ double norm = sqrt(nx*nx+ny*ny);
+ nx = nx / norm;
+ ny = ny / norm;
+
+ // Computation of the integration points (parametric) and weights.
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
+
+ // Computation of the sum (of integrand evaluated at integration point * weight)
+ double gSum = 0;
+ for(size_t k = 0; k < weights.size(); k++)
+ {
+ // Parametric coordinate
+ double integration_variable = localCoord[3*k];
+ // Global coordinates
+ double hX = element.midX + element.deltaX*integration_variable/2.0;
+ double hY = element.midY + element.deltaY*integration_variable/2.0;
+
+ // Distance between mid point of element 'i' and integration point.
+ double r = sqrt(pow(x - hX, 2.0) + pow(y - hY, 2.0));
+ gSum += (-((1.0/2.0*pow(r, 2) - pow((xi - hX),2))/(pow(r, 4)))*nx + ((xi - hX)*(yi - hY)*ny)/(pow(r, 4))) * weights[k];
+ }
+
+ double gradHx = gSum * element.length / (2*M_PI);
+
+ return gradHx;
+}
+double computeNumericalGradHy(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType)
+{
+ double xi = x;
+ double yi = y;
+ double xj = element.x1;
+ double xjp1 = element.x2;
+ double yj = element.y1;
+ double yjp1 = element.y2;
+
+ double dx = xjp1-xj;
+ double dy = yjp1-yj;
+ double nx = dy;
+ double ny = -dx;
+ double norm = sqrt(nx*nx+ny*ny);
+ nx = nx / norm;
+ ny = ny / norm;
+
+ // Computation of the integration points (parametric) and weights.
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
+
+ // Computation of the sum (of integrand evaluated at integration point * weight)
+ double gSum = 0;
+ for(size_t k = 0; k < weights.size(); k++)
+ {
+ // Parametric coordinate
+ double integration_variable = localCoord[3*k];
+ // Global coordinates
+ double hX = element.midX + element.deltaX*integration_variable/2.0;
+ double hY = element.midY + element.deltaY*integration_variable/2.0;
+
+ // Distance between mid point of element 'i' and integration point.
+ double r = sqrt(pow(xi - hX, 2.0) + pow(yi - hY, 2.0));
+ gSum += (-((1.0/2.0*pow(r, 2) - pow((yi - hY),2))/(pow(r, 4)))*ny + ((xi - hX)*(yi - hY)*nx)/(pow(r, 4))) * weights[k];
+ }
+
+ double gradHy = gSum * element.length / (2*M_PI);
+
+ return gradHy;
+}
+double computeNumericalGradGx(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType)
+{
+ // Computation of the integration points (parametric) and weights.
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
+
+ // Computation of the sum (of integrand evaluated at integration point * weight)
+ double gSum = 0;
+ for(size_t k = 0; k < weights.size(); k++)
+ {
+ // Parametric coordinate
+ double integration_variable = localCoord[3*k];
+ // Global coordinates
+ double gX = element.midX + element.deltaX*integration_variable/2.0;
+ double gY = element.midY + element.deltaY*integration_variable/2.0;
+
+ // Distance between mid point of element 'i' and integration point.
+ double r = sqrt(pow(x - gX, 2.0) + pow(y - gY, 2.0));
+ gSum += ((x-gX)/(pow(r, 2))) * weights[k];
+ }
+
+ double gradGx = gSum * element.length / (4*M_PI);
+
+ return gradGx;
+}
+double computeNumericalGradGy(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType)
+{
+ // Computation of the integration points (parametric) and weights.
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
+
+ // Computation of the sum (of integrand evaluated at integration point * weight)
+ double gSum = 0;
+ for(size_t k = 0; k < weights.size(); k++)
+ {
+ // Parametric coordinate
+ double integration_variable = localCoord[3*k];
+ // Global coordinates
+ double gX = element.midX + element.deltaX*integration_variable/2.0;
+ double gY = element.midY + element.deltaY*integration_variable/2.0;
+
+ // Distance between mid point of element 'i' and integration point.
+ double r = sqrt(pow(x - gX, 2.0) + pow(y - gY, 2.0));
+ gSum += ((y-gY)/(pow(r, 2))) * weights[k];
+ }
+
+ double gradGy = gSum * element.length / (4*M_PI);
+
+ return gradGy;
+}
\ No newline at end of file
diff --git a/srcs/BEM/functionsBEM.cpp b/srcs/BEM/functionsBEM.cpp
deleted file mode 100644
index 0203b1ac7538aaccecaeaa475f4da3c25ed5c72a..0000000000000000000000000000000000000000
--- a/srcs/BEM/functionsBEM.cpp
+++ /dev/null
@@ -1,1230 +0,0 @@
-#include <gmsh.h>
-#include <iostream>
-#include <map>
-#include <sstream>
-#include <Eigen/Dense>
-#include <cmath>
-#include <iomanip>
-#include <stdio.h>
-#include <chrono>
-#include "functionsBEM.hpp"
-#include <omp.h>
-#include <algorithm>
-#include <list>
-
-/*-----------------todo déjà fait mais il faudra peut être en refaire certains par après donc je les garde------------------*/
-//todo: séparer en différentes fonctions
-//todo: comparer au livre BEM
-//todo: voir pour gérer des trous dans les surfaces
-//todo: combiner avec la FEM
-//todo: faire plein de tests différent (dans les .geo)
-//todo: afficher et calculer la force electrostatique et champs electriques
-//todo: on a besoin de la pression electrostatique à la frontière mais on peut essayer de l'afficher partout dans le domaine => representation formula pour le domain ou par elementnodedata et gradient
-//todo: fusionner avec PY
-//todo: potentiel2d fonction
-//todo: gérer permuttation matrice plus efficacement (en vrai mnt, osef)
-//todo: hybrid mesh + éléments 1d 3 noeuds
-//todo: mettre en parallèle
-//todo: elementtype pour optimiser elementnodedata
-//todo: commenter le code
-//todo: faire une fonction pour chrono afin de mesurer le temps de toutes les différentes parties du code
-//todo: convergence => utiliser une erreur L2
-//todo: gérer des matériaux composites
-//todo: implémenter solver dynamique
-//todo: demander pour le critical d'openmp
-
-/*-----------------todo qu'il reste vraiment à faire------------------*/
-//todo: afficher les déplacements, le potentiel en fonction de la distance dans l'élément (une coupe 1D dans le domaine)
-//todo: representation formula pour le champ electrique
-//todo: créer d'autres .geo avec des géométries plus complexes
-
-//astuce: pour reverse le sens de parcours d'une courbe et inverser les tangentes dans le .geo, utiliser reversemesh
-
-//todo: openmp
-//todo: séparer analytique et numérique
-//todo: visu electric field, combiner avec les autres
-//todo: refaire le vecteur data pour elementdata (plus d'elementcount)
-
-void element_node_data(std::vector<std::vector<std::vector<double>>> &data, std::vector<std::vector<size_t>> &elementTags2d_domain_BEM, std::vector<int> &elementTypes2d_domain_BEM, std::vector<std::vector<size_t>> &nodeTags2d_domain_BEM, std::vector<elementStruct> &elementVector, std::map<int, double> &nodeMapDomain, bool &show_element_node_data, int &entities_tags_for_physical_group_domain_BEM, std::map<std::size_t, std::pair<bool,bool>> &nodeOnBoundary)
-{
- //The first part of the algorithm consists in recovering the tag of the 2d elements of the domain which have at least one node in common with the boundary
- //We initialise the first dimension of the vector that will contain all 2d elements in common with the boundary with the number of element types present in the domain.
- std::vector<std::vector<size_t>> elementTags2d_boundary(elementTypes2d_domain_BEM.size());
-
- //this loop goes through all the nodes of the boundary (node 1 of each element). We then use getElementsByCoordinates to obtain all the elements in common with the node considered. We only take the 2d elements.
- for(size_t i = 0; i < elementVector.size(); i++)
- {
- double posx = elementVector[i].x1;
- double posy = elementVector[i].y1;
- double posz = 0;
- std::vector<size_t> elementTags;
- gmsh::model::mesh::getElementsByCoordinates(posx,posy,posz,elementTags,2);
-
- //we loop on the 2d elements which are in common with the considered node of the boundary. Depending on their type, they are added to the vector comprising all the 2d elements common with the boundary.
- for(size_t j = 0; j < elementTags.size(); j++)
- {
- int elementType;
- std::vector<size_t> nodeTags;
- int dim;
- int tag;
- std::vector<int> physicalTags;
- gmsh::model::mesh::getElement(elementTags[j],elementType,nodeTags,dim,tag);
- gmsh::model::getPhysicalGroupsForEntity(dim,tag,physicalTags);
-
- for(size_t k = 0; k < elementTypes2d_domain_BEM.size(); k++)
- {
- for(size_t l = 0; l < physicalTags.size(); l++)
- {
- if(elementType == elementTypes2d_domain_BEM[k] && physicalTags[l] == entities_tags_for_physical_group_domain_BEM)
- {
- elementTags2d_boundary[k].push_back(elementTags[j]);
- }
- }
- }
- }
- }
-
- //This loop allows duplicate tags to be removed. If an element has several nodes with the border, it has been counted several times before.
- for(size_t i = 0; i < elementTypes2d_domain_BEM.size(); i++)
- {
- sort( elementTags2d_boundary[i].begin(), elementTags2d_boundary[i].end() );
- elementTags2d_boundary[i].erase( unique( elementTags2d_boundary[i].begin(), elementTags2d_boundary[i].end() ), elementTags2d_boundary[i].end() );
- }
-
- //This map may be replaced in the future. It allows to retain the position of the elements of the elementTag2d loop and linked their position to the loop of the elementTag2d_boundary (the tags of the elements in common with the boundary). This map is the link between the array of tags of the 2d elements of the boundary and that of the tags of all the 2d elements.
- std::map<int, int> tagposition;
-
- //The first dimension of data corresponds to the number of element types.
- data.resize(elementTags2d_domain_BEM.size());
- for(size_t i = 0; i < elementTags2d_domain_BEM.size(); i++)
- {
- //The second dimension of data corresponds to the number of element tags for on element type.
- data[i].resize(elementTags2d_domain_BEM[i].size());
- std::string elementName;
- int dim;
- int order;
- int numNodes;
- std::vector<double> localNodeCoord;
- int numPrimaryNodes;
- gmsh::model::mesh::getElementProperties(elementTypes2d_domain_BEM[i], elementName, dim, order, numNodes, localNodeCoord, numPrimaryNodes);
-
- //getElementProperties was used to get the number of nodes for a specific element
- int nb_nodes_element2d = numNodes;
- for(size_t j = 0; j < elementTags2d_domain_BEM[i].size(); j++)
- {
- //The third dimension of data corresponds to the number of nodes(third: k loop) for an element tags(second: j loop) for an element type(first: i loop).
- data[i][j].resize(nb_nodes_element2d,HUGE_VAL);
- tagposition[elementTags2d_domain_BEM[i][j]] = j;
- }
- }
-
- //The second part of the algorithm will consist of assigning the correct potential value to the nodes of the elements in common with the boundary.
- //The first loop is on the elementType in 2d in common with the boundary.
- for(size_t i = 0; i < elementTags2d_boundary.size(); i++)
- {
- std::string elementName;
- int dim;
- int order;
- int numNodes;
- std::vector<double> localNodeCoord;
- int numPrimaryNodes;
- gmsh::model::mesh::getElementProperties(elementTypes2d_domain_BEM[i], elementName, dim, order, numNodes, localNodeCoord, numPrimaryNodes);
- //getElementProperties was used to get the number of nodes for a specific element
- int nb_nodes_element2d = numNodes;
-
- //The second loop is on the elementTag in 2d in common with the boundary.
- for(size_t j = 0; j < elementTags2d_boundary[i].size(); j++)
- {
- std::vector<int> elements_commun_tag(0,0);
- std::vector<int> elements_commun_indice(0,0);
- std::vector<int> nodes_commun_tag(0,0);
- int elementType;
- std::vector<size_t> nodeTags_element2d;
- int dim;
- int tag;
- gmsh::model::mesh::getElement(elementTags2d_boundary[i][j],elementType,nodeTags_element2d,dim,tag);
- //The third loop is on the nodes of a specific elementTag in 2d in common with the boundary.
- for(int k = 0; k < nb_nodes_element2d; k++)
- {
- int loop_element_2d_nodes_tag = nodeTags_element2d[k];
- //For each node of the 2d elements, we loop over all the 1d elements of the boundary and look for both the node in common and the elements in common. Indeed, on the boundary, a specific node belongs to two 1d elements.
- for(size_t l = 0; l < elementVector.size(); l++)
- {
- //condition vérifié si le noeud de l’élément 2d est égal au noeud d’un élément 1d de la frontière.
- if(elementVector[l].nodeTags[0] == (size_t) loop_element_2d_nodes_tag || elementVector[l].nodeTags[1] == (size_t) loop_element_2d_nodes_tag)
- {
- //this loop allows the node to be added only once if it is present several times. As the 1d elements of the boundary are traversed, 1 node encounters 2 elements and can therefore be counted multiple times.
- int compteur_node = 0;
- for(size_t m = 0; m < nodes_commun_tag.size(); m++)
- {
- if(nodes_commun_tag[m] == loop_element_2d_nodes_tag)
- {
- compteur_node++;
- }
- }
- if(compteur_node == 0)
- {
- nodes_commun_tag.push_back(loop_element_2d_nodes_tag);
- }
- //The same procedure is applied for the elements. You only want to count them once if they are in common
- int compteur_element = 0;
- for(size_t m = 0; m < elements_commun_tag.size(); m++)
- {
- if(elements_commun_tag[m] == elementVector[l].tag)
- {
- compteur_element++;
- }
- }
- if(compteur_element == 0)
- {
- //Both the tag of the element and its position in the loop on the elements of the boundary are retained.
- elements_commun_tag.push_back(elementVector[l].tag);
- elements_commun_indice.push_back(l);
- }
- }
- }
- }
- //For each 2d element in common with the boundary, displays the tag of the nodes and the tag of the common 1d elements of the boundary.
- if(show_element_node_data)
- {
- std::cout << "noeuds commun: "<< nodes_commun_tag.size() << ": ";
- for(size_t m = 0; m < nodes_commun_tag.size(); m++)
- {
- std::cout << nodes_commun_tag[m] << " ";
- }
- std::cout << std::endl;
- std::cout << "elements commun: " << elements_commun_tag.size() << ": ";
- for(size_t m = 0; m < elements_commun_tag.size(); m++)
- {
- std::cout << elements_commun_tag[m] << " ";
- }
- std::cout << std::endl;
- std::cout << std::endl;
- }
- //Several geometric configurations are to be considered. Here, the 2d element has 1 node in common with the boundary and must therefore touch 2 of its elements.
- if(nodes_commun_tag.size() == 1)
- {
- if(elements_commun_tag.size() == 2)
- {
- //We loop over the nodes of the 2d element, find the one that is common to the boundary and assign it the average value of the potential of the 2 elements it touches. Indeed the potential is discontinuous from one element of the boundary to the other.
- for(int k = 0; k < nb_nodes_element2d; k++)
- {
- int loop_element_2d_nodes_tag = nodeTags_element2d[k];
- if(loop_element_2d_nodes_tag == nodes_commun_tag[0])
- {
- double phi = 0;
- phi = elementVector[elements_commun_indice[0]].bcValue[0];
- phi = phi + elementVector[elements_commun_indice[1]].bcValue[0];
- phi = phi/2.0;
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = phi;
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 0)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from common nodes: " << nodes_commun_tag.size() << " and common elements: " << elements_commun_tag.size() << std::endl;
- }
- else
- {
- nodeOnBoundary[loop_element_2d_nodes_tag].second = 1;
- }
- }
- }
- }
- //displays an error in the case of geometry where 1 node is in contact with more than 2 elements.
- else
- {
- disp_element_node_data(nodes_commun_tag, elements_commun_tag);
- }
- }
- //The principle is the same as before. For the case where there are 2 nodes in common, the 2d element can touch either 3 or 4 1d elements of the boundary. In the first case, it has one of its sides in common with the boundary, the other case can be found in the corners of the domain (two of its corners touch the boundary in 1 node and meet 2 elements like case 1).
- if(nodes_commun_tag.size() == 2)
- {
- if(elements_commun_tag.size() == 3)
- {
- for(int k = 0; k < nb_nodes_element2d; k++)
- {
- int loop_element_2d_nodes_tag = nodeTags_element2d[k];
- for(size_t l = 0; l < elements_commun_tag.size(); l++)
- {
- //Among the 3 elements of the boundary in common, we look for the one which has 2 nodes in common with the 2d element. We look for the side and therefore the 1d element that they have in common.
- int compteur_nodes = 0;
- for(size_t m = 0; m < nodes_commun_tag.size(); m++)
- {
- if(elementVector[elements_commun_indice[l]].nodeTags[0] == (size_t) nodes_commun_tag[m]||elementVector[elements_commun_indice[l]].nodeTags[1] == (size_t) nodes_commun_tag[m])
- {
- compteur_nodes++;
- }
- }
- if(compteur_nodes == 2)
- {
- //Once we have found the element in common, the node of the 2d element on which we are looping must be the right one
- if((size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].nodeTags[0]|| (size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].nodeTags[1])
- {
- //The constant value of the boundary element in common is assigned to the 2 nodes of the 2d element that touches the boundary.
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = elementVector[elements_commun_indice[l]].bcValue[0];
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 0)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from common nodes: " << nodes_commun_tag.size() << " and common elements: " << elements_commun_tag.size() << std::endl;
- }
- else
- {
- nodeOnBoundary[loop_element_2d_nodes_tag].second = 1;
- }
- }
- //In the case of second order elements, the value of the common 1d boundary element is assigned to the middle node.
- if(order == 2)
- {
- if((size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].midNodeTags)
- {
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = elementVector[elements_commun_indice[l]].bcValue[0];
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 0)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from common nodes: " << nodes_commun_tag.size() << " and common elements: " << elements_commun_tag.size() << std::endl;
- }
- else
- {
- nodeOnBoundary[loop_element_2d_nodes_tag].second = 1;
- }
- }
- }
- }
- }
- }
- }
- else if(elements_commun_tag.size() == 4)
- {
- for(int k = 0; k < nb_nodes_element2d; k++)
- {
- int loop_element_2d_nodes_tag = nodeTags_element2d[k];
- for(size_t l = 0; l < elements_commun_tag.size(); l++)
- {
- //This geometry behaves almost like the first case (1 node, 2 elements) except that it is true for 2 different nodes of the 2d element. We loop over the elements of the boundary. If one of its nodes is common with the boundary, it is assigned the value of the potential of this first element. When the second one comes, it is averaged with the potential already present on the node.
- if((size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].nodeTags[0]|| (size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].nodeTags[1])
- {
- if(data[i][tagposition[elementTags2d_boundary[i][j]]][k] == HUGE_VAL)
- {
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = elementVector[elements_commun_indice[l]].bcValue[0];
- }
- else
- {
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = (data[i][tagposition[elementTags2d_boundary[i][j]]][k] + elementVector[elements_commun_indice[l]].bcValue[0])/2.0;
- }
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 0)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from common nodes: " << nodes_commun_tag.size() << " and common elements: " << elements_commun_tag.size() << std::endl;
- }
- else
- {
- nodeOnBoundary[loop_element_2d_nodes_tag].second = 1;
- }
- }
- }
- }
- }
- else
- {
- disp_element_node_data(nodes_commun_tag, elements_commun_tag);
- }
- }
- //The last case is that of 3 common nodes which therefore affect 4 elements. This case is mainly found in the corners. The average potential is attributed to the node that touches 2 elements (in the corner of a square for example). The other 2 nodes must simply match both the side in common and the node in common. The principle of assigning the potential value of the node is the superposition of the previous methods.
- if(nodes_commun_tag.size() == 3)
- {
- if(elements_commun_tag.size() == 4)
- {
- for(int k = 0; k < nb_nodes_element2d; k++)
- {
- int loop_element_2d_nodes_tag = nodeTags_element2d[k];
- for(size_t l = 0; l < elements_commun_tag.size(); l++)
- {
- int compteur_nodes = 0;
- for(size_t m = 0; m < nodes_commun_tag.size(); m++)
- {
- if(elementVector[elements_commun_indice[l]].nodeTags[0] == (size_t) nodes_commun_tag[m]||elementVector[elements_commun_indice[l]].nodeTags[1] == (size_t) nodes_commun_tag[m])
- {
- compteur_nodes++;
- }
- }
- if(compteur_nodes == 2)
- {
- if((size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].nodeTags[0]|| (size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].nodeTags[1])
- {
- if(data[i][tagposition[elementTags2d_boundary[i][j]]][k] == HUGE_VAL)
- {
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = elementVector[elements_commun_indice[l]].bcValue[0];
- }
- else
- {
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = (data[i][tagposition[elementTags2d_boundary[i][j]]][k] + elementVector[elements_commun_indice[l]].bcValue[0])/2.0;
- }
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 0)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from common nodes: " << nodes_commun_tag.size() << " and common elements: " << elements_commun_tag.size() << std::endl;
- }
- else
- {
- nodeOnBoundary[loop_element_2d_nodes_tag].second = 1;
- }
- }
- //In the case of second order elements, the value of the common 1d boundary element is assigned to the middle node.
- if(order == 2)
- {
- if((size_t) loop_element_2d_nodes_tag == elementVector[elements_commun_indice[l]].midNodeTags)
- {
- data[i][tagposition[elementTags2d_boundary[i][j]]][k] = elementVector[elements_commun_indice[l]].bcValue[0];
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 0)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from common nodes: " << nodes_commun_tag.size() << " and common elements: " << elements_commun_tag.size() << std::endl;
- }
- else
- {
- nodeOnBoundary[loop_element_2d_nodes_tag].second = 1;
- }
- }
- }
- }
- }
- }
- }
- else
- {
- disp_element_node_data(nodes_commun_tag, elements_commun_tag);
- }
- }
- //We clear the vectors that we push at each loop
- elements_commun_tag.clear();
- elements_commun_indice.clear();
- nodes_commun_tag.clear();
- }
- }
-
-
- //The last step is to loop over all the 2d elements in the domain. If the HUGE_VAL value is still present, this means that the element node does not belong to the boundary (since the value of the potential has not been modified) and is therefore part of the domain. We then use the map which associates the value of the potential with the tag of the node inside the domain.
- for(size_t i = 0; i < elementTags2d_domain_BEM.size(); i++)
- {
- std::string elementName;
- int dim;
- int order;
- int numNodes;
- std::vector<double> localNodeCoord;
- int numPrimaryNodes;
- gmsh::model::mesh::getElementProperties(elementTypes2d_domain_BEM[i], elementName, dim, order, numNodes, localNodeCoord, numPrimaryNodes);
- int nb_nodes_element2d = numNodes;
- for(size_t j = 0; j < elementTags2d_domain_BEM[i].size(); j++)
- {
- int elementType;
- std::vector<size_t> nodeTags_element2d;
- int dim;
- int tag;
- gmsh::model::mesh::getElement(elementTags2d_domain_BEM[i][j],elementType,nodeTags_element2d,dim,tag);
- for(int k = 0; k < nb_nodes_element2d; k++)
- {
- int loop_element_2d_nodes_tag = nodeTags_element2d[k];
- //je laisse la condition hugeval pour l'instant car elle est plus "gentille". Le code se plantera moins facilement mais il y aura quand même un message dans la console si on a des noeuds qui ont un potentiel anormal
- if(data[i][j][k] == HUGE_VAL)
- {
- if(nodeMapDomain.count(loop_element_2d_nodes_tag))
- {
- data[i][j][k] = nodeMapDomain[loop_element_2d_nodes_tag];
- }
- else
- {
- std::cout << "error: node not found in nodeMapDomain code BEM" << std::endl;
- }
- if(nodeOnBoundary[loop_element_2d_nodes_tag].first == 1 || nodeOnBoundary[loop_element_2d_nodes_tag].second == 1)
- {
- std::cerr << "error node: " << loop_element_2d_nodes_tag << " is considered on the boundary and in the BEM domain" << " output from the use of nodemapdomain" << std::endl;
- }
- }
- }
- }
- }
-}
-
-void dispNodeCoordFun(const std::vector<size_t> &allNodeTags, const std::vector<double> &allCoord)
-{
- for(size_t i = 0; i < allNodeTags.size(); i++)
- {
- std::cout << "Node '" << allNodeTags[i] << "': (" << allCoord[3*i] << "; " << allCoord[3*i + 1] << "; " << allCoord[3*i + 2] << ")\n";
- }
- std::cout << std::endl;
-}
-
-void dispMatricesFun(const Eigen::MatrixXd &A, const Eigen::MatrixXd &B, const Eigen::MatrixXd &c, const Eigen::MatrixXd &b, const Eigen::MatrixXd &x)
-{
- std::cout << "A matrix:\n";
- std::cout << A << std::endl;
- std::cout << "\nB matrix:\n";
- std::cout << B << std::endl;
- std::cout << "\nc vector:\n";
- std::cout << c << std::endl;
- std::cout << "\nb vector:\n";
- std::cout << b << std::endl;
- std::cout << "\nx vector:\n";
- std::cout << x << std::endl;
-}
-
-void dispValuesFun(const std::vector<elementStruct> &elementVector)
-{
- for(size_t i = 0; i < elementVector.size(); i++)
- {
- std::cout << "Element '" << elementVector[i].tag << "': from node " << elementVector[i].nodeTags[0] << " to node " << elementVector[i].nodeTags[1] << " (" << elementVector[i].bcType << "):\tDirichlet: " << elementVector[i].bcValue[0] << "\tNeumann: " << elementVector[i].bcValue[1] << "\n";
- }
-}
-
-void getBC(std::map<size_t, int> &allNodeIndex, const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy, std::vector<size_t> &NodeTags2d_domain, std::vector<double> &NodeCoord2d_domain, std::vector<int> &elementTypes2d_domain_BEM, std::vector<std::vector<size_t>> &elementTags2d_domain_BEM, std::vector<std::vector<size_t>> &elnodeTags2d_domain_BEM, int &entities_tags_for_physical_group_domain_BEM, std::vector<int> &entities_tags_2d_domain_BEM, std::map<int, std::pair<double, double>> &boundaryDisplacementMap, std::map<std::size_t, std::pair<bool,bool>> &nodeOnBoundary, const int domainTag)
-{
- // pour moi nodeTags2d_domain inutile, ainsi que nodeCoords2d_domain
- // Gets the dim and the tags of the physical groups
- gmsh::vectorpair physicalGroupsDimTags;
- gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
-
- // Gets the boundary conditions from onelab and create a map: "name of the physical group" -> ("dirichlet"/"neumann", value)
- std::map<std::string, std::pair<std::string, double>> bcMap;
-
- std::vector<std::string> bcKeys;
- gmsh::onelab::getNames(bcKeys, "Boundary Conditions.+");
-
- std::string domainName = "BEM_domain_" + std::to_string(domainTag);
- for(auto &key : bcKeys)
- {
- std::vector<double> values;
- gmsh::onelab::getNumber(key, values);
-
- std::stringstream ss(key);
- std::vector<std::string> words;
- std::string word;
- while(getline(ss, word, '/'))
- words.push_back(word);
-
- // focus only on dirichlet or neumann boundary conditions -> associated to BEM // CHECK BON DOMAIN
- if(words.size() == 4 && (words[3].compare("neumann") == 0 || words[3].compare("dirichlet") == 0) && words[2].compare(domainName) == 0)
- {
- bcMap[words[1]] = {words[3], values[0]}; // doit etre particularisée
- }
- }
-
- // Initially, these are 0 and are incremented each time a dirichlet/neumann-element is encountered.
- for (size_t i = 0; i < physicalGroupsDimTags.size(); i++)
- {
- int physicalGroupsDim = physicalGroupsDimTags[i].first;
- int physicalGroupsTag = physicalGroupsDimTags[i].second;
- std::string name;
- gmsh::model::getPhysicalName(physicalGroupsDim, physicalGroupsTag, name);
-
- if(name == domainName)
- {
- gmsh::model::mesh::getNodesForPhysicalGroup(physicalGroupsDim,physicalGroupsTag,NodeTags2d_domain,NodeCoord2d_domain);
-
- entities_tags_for_physical_group_domain_BEM = physicalGroupsTag;
- gmsh::model::getEntitiesForPhysicalGroup(physicalGroupsDim, physicalGroupsTag, entities_tags_2d_domain_BEM);
-
- for(size_t j=0; j< entities_tags_2d_domain_BEM.size(); j++)
- {
- std::vector<int> elementTypes;
- std::vector<std::vector<size_t> > elementTags;
- std::vector<std::vector<size_t> > nodeTags;
- gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, physicalGroupsDim, entities_tags_2d_domain_BEM[j]);
- elementTypes2d_domain_BEM.insert(elementTypes2d_domain_BEM.end(), elementTypes.begin(), elementTypes.end());
- elementTags2d_domain_BEM.insert(elementTags2d_domain_BEM.end(), elementTags.begin(), elementTags.end());
- }
- //remove_no_domain_nodes(NodeTags2d_domain, elementVector); //ça fait bugger
- }
-
- //s'il ne trouve pas name dans la map, find retourne map.end() => condition pour continue
- if(bcMap.find(name) == bcMap.end())
- {
- continue;
- }
-
- // Retrieves the type and the value of the boundary condition of this physical group.
- // /!\ Should check that this physical groups was well defined in the list of "Boundary Condition" !
- std::string bcType = bcMap[name].first;
- double bcValue = bcMap[name].second;
-
- // Prints the physical groups name, dim, tag, bcType and bcValue (if dispHierarchy != 0).
- if(dispHierarchy)
- std::cout << "Physical group '" << name << "' (" << physicalGroupsDim << "D) with tag '" << physicalGroupsTag << "' has a '" << bcType << "' type and the value " << bcValue << "\n";
-
- std::vector<int> entitiesTags;
- gmsh::model::getEntitiesForPhysicalGroup(physicalGroupsDim, physicalGroupsTag, entitiesTags);
-
- // devrait passer seulement sur les physical groups qui sont à la frontière de BEM_domain_i + Neumann ou Dirichlet
- for (size_t k = 0; k < entitiesTags.size(); k++)
- {
- // Print the tag of the entities of each physical group (if dispHierarchy != 0).
- if(dispHierarchy)
- std::cout << "\t contains entity '" << entitiesTags[k] << "' ";
-
-
- std::vector<size_t> elementTags1;
- std::vector<size_t> nodeTags1;
- std::vector<size_t> elementTags8;
- std::vector<size_t> nodeTags8;
-
- gmsh::model::mesh::getElementsByType(1, elementTags1, nodeTags1, entitiesTags[k]);
- gmsh::model::mesh::getElementsByType(8, elementTags8, nodeTags8, entitiesTags[k]);
-
- // Print the number of element of the given entity (if dispHierarchy != 0).
- if(dispHierarchy)
- std::cout << "which has " << elementTags1.size() + elementTags8.size() << " elements:\n";
-
- for(size_t el = 0; el < elementTags1.size(); el++)
- {
- int elTag = elementTags1[el];
-
- elementStruct element;
- element.tag = elTag;
- element.nodeTags[0] = nodeTags1[2*el];
- element.nodeTags[1] = nodeTags1[2*el + 1];
- element.midNodeTags = 0;
-
- std::pair<double, double> nodalDisp1 = {0.0, 0.0}, nodalDisp2 = {0.0, 0.0};
- if(boundaryDisplacementMap.count(element.nodeTags[0]))
- nodalDisp1 = boundaryDisplacementMap[element.nodeTags[0]];
- if(boundaryDisplacementMap.count(element.nodeTags[1]))
- nodalDisp2 = boundaryDisplacementMap[element.nodeTags[1]];
- computeGeometricParam(element, allNodeIndex, allCoord, nodalDisp1, nodalDisp2);
-
- if(bcType == "dirichlet")
- {
- element.bcValue[0] = bcValue;
- elementVector.insert(elementVector.begin(), element);
- dirichletCount++;
- }
- else if(bcType == "neumann")
- {
- element.bcValue[1] = bcValue;
- elementVector.push_back(element);
- }
- // Print the tag of a given element and the tags of the nodes that composes it as well as their coordinates (if dispHierarchy != 0).
- if(dispHierarchy)
- {
- std::cout << "\t\t-Element '" << elTag << "' (" << element.bcType << "):\n";
- std::cout << "\t\t\t-> Node '" << element.nodeTags[0] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[0]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[0]] + 1] << ")\n";
- std::cout << "\t\t\t-> Node '" << element.nodeTags[1] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[1]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[1]] + 1] << ")\n";
- }
- }
- for(size_t el = 0; el < elementTags8.size(); el++)
- {
- int elTag = elementTags8[el];
-
- elementStruct element;
- element.tag = elTag;
- element.nodeTags[0] = nodeTags8[3*el];
- element.nodeTags[1] = nodeTags8[3*el + 1];
- element.midNodeTags = nodeTags8[3*el + 2];
-
- std::pair<double, double> nodalDisp1 = {0.0, 0.0}, nodalDisp2 = {0.0, 0.0};
- if(boundaryDisplacementMap.count(element.nodeTags[0]))
- nodalDisp1 = boundaryDisplacementMap[element.nodeTags[0]];
- if(boundaryDisplacementMap.count(element.nodeTags[1]))
- nodalDisp2 = boundaryDisplacementMap[element.nodeTags[1]];
- computeGeometricParam(element, allNodeIndex, allCoord, nodalDisp1, nodalDisp2);
-
- if(bcType == "dirichlet")
- {
- element.bcValue[0] = bcValue;
- elementVector.insert(elementVector.begin(), element);
- dirichletCount++;
- }
- else if(bcType == "neumann")
- {
- element.bcValue[1] = bcValue;
- elementVector.push_back(element);
- }
-
- // Print the tag of a given element and the tags of the nodes that composes it as well as their coordinates (if dispHierarchy != 0).
- if(dispHierarchy)
- {
- std::cout << "\t\t-Element '" << elTag << "' (" << element.bcType << "):\n";
- std::cout << "\t\t\t-> Node '" << element.nodeTags[0] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[0]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[0]] + 1] << ")\n";
- std::cout << "\t\t\t-> Node '" << element.nodeTags[1] << "' (" << allCoord[3*allNodeIndex[element.nodeTags[1]]] << "; " << allCoord[3*allNodeIndex[element.nodeTags[1]] + 1] << ")\n";
- }
- }
- if(dispHierarchy)
- std::cout << "\n";
- }
- }
- for(std::size_t i = 0; i < NodeTags2d_domain.size(); i++)
- {
- nodeOnBoundary[NodeTags2d_domain[i]].first = 0;
- nodeOnBoundary[NodeTags2d_domain[i]].second = 0;
- }
- for(std::size_t i = 0; i < elementVector.size(); i++)
- {
- nodeOnBoundary[elementVector[i].nodeTags[0]].first = 1;
- nodeOnBoundary[elementVector[i].nodeTags[1]].first = 1;
- if(elementVector[i].midNodeTags != 0)
- {
- nodeOnBoundary[elementVector[i].midNodeTags].first = 1;
- }
- }
-}
-
-void computeGeometricParam(elementStruct &element, std::map<std::size_t, int> &allNodeIndex, const std::vector<double> &allCoord, std::pair<double, double> &nodalDisp1, std::pair<double, double> &nodalDisp2)
-{
- int i1Tag = element.nodeTags[0];
- int i2Tag = element.nodeTags[1];
-
- int i1Index = allNodeIndex[i1Tag];
- int i2Index = allNodeIndex[i2Tag];
-
- // Get the coordinates of the nodes from their tag.
- element.x1 = allCoord[3*i1Index] + nodalDisp1.first;
- element.y1 = allCoord[3*i1Index + 1] + nodalDisp1.second;
- element.x2 = allCoord[3*i2Index] + nodalDisp2.first;
- element.y2 = allCoord[3*i2Index + 1] + nodalDisp2.second;
-
- element.midX = (element.x1 + element.x2) / 2;
- element.midY = (element.y1 + element.y2) / 2;
- element.deltaX = element.x2 - element.x1;
- element.deltaY = element.y2 - element.y1;
- element.length = sqrt(pow(element.deltaX, 2.0) + pow(element.deltaY, 2.0));
-}
-
-
-double computeH(const double &x, const double &y, const elementStruct &element)
-{
- double deltaX1 = element.x1 - x;
- double deltaX2 = element.x2 - x;
- double deltaY1 = element.y1 - y;
- double deltaY2 = element.y2 - y;
- double cosAlpha = (deltaX1*deltaX2 + deltaY1*deltaY2) / (sqrt((pow(deltaX1,2.0) + pow(deltaY1,2.0)) * (pow(deltaX2,2.0) + pow(deltaY2,2.0))));
-
- if(cosAlpha < -1.0)
- {
- cosAlpha = -1.0;
- }
- else if(cosAlpha > 1.0)
- {
- cosAlpha = 1.0;
- }
-
- double alpha = acos(cosAlpha);
- double crossProduct = deltaX1*deltaY2 - deltaY1*deltaX2;
-
- if(crossProduct < 0)
- {
- alpha = -alpha;
- }
-
- return alpha / (2*M_PI);
-}
-
-
-double computeG(const double &x, const double &y, const elementStruct &element, const std::string &integrationType)
-{
- // Computation of the integration points (parametric) and weights.
- std::vector<double> localCoord;
- std::vector<double> weights;
- gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
-
- // Computation of the sum (of integrand evaluated at integration point * weight)
- double gSum = 0;
- for(size_t k = 0; k < weights.size(); k++)
- {
- // Parametric coordinate
- double xi = localCoord[3*k];
- // Global coordinates
- double gX = element.midX + element.deltaX*xi/2;
- double gY = element.midY + element.deltaY*xi/2;
-
- // Distance between mid point of element 'i' and integration point.
- double r = sqrt(pow(x - gX, 2.0) + pow(y - gY, 2.0));
- gSum += log(r) * weights[k];
- }
-
- double g = gSum * element.length / (4*M_PI);
-
- return g;
-}
-
-void displayData(const std::string viewName, const std::string modelName, const std::string fileName, const std::vector<size_t> &tags, const std::vector<std::vector<double>> &data)
-{
- int viewTag = gmsh::view::add(viewName);
- gmsh::view::addModelData(viewTag, 0, modelName, "ElementData", tags, data);
- gmsh::view::write(viewTag, fileName);
-
- //gmsh::option::setNumber("View[0].IntervalsType", 3);
- // option::setNumber("View[0].AdaptVisualizationGrid", 1); // Problem with this command !
- //gmsh::option::setNumber("View[0].MaxRecursionLevel", 3);
- //gmsh::option::setNumber("View[0].TargetError", -0.0001);
-}
-
-void compute_internal_phi(std::vector<double> &allCoord, std::vector<size_t> &NodeTags2d_domain, std::map<size_t, int> &allNodeIndex, std::vector<elementStruct> &elementVector, std::map<int,double> &nodeMapDomain, const std::string integrationType, std::map<std::size_t, std::pair<bool,bool>> &nodeOnBoundary, bool &solution_type)
-{
- double x = 0;
- double y = 0;
- double phi = 0;
- double h = 0;
- double g = 0;
- std::vector<double> vectorForNodeMap(NodeTags2d_domain.size(),0);
-
-#pragma omp parallel for schedule(static) private(x,y,phi,h,g) shared(NodeTags2d_domain,allCoord,elementVector,nodeMapDomain,nodeOnBoundary,solution_type)
- for(std::size_t i = 0; i < NodeTags2d_domain.size(); i++)
- {
- if(nodeOnBoundary[NodeTags2d_domain[i]].first == 0)
- {
- // Coordinates of the node
- x = allCoord[3*allNodeIndex[NodeTags2d_domain[i]]];
- y = allCoord[3*allNodeIndex[NodeTags2d_domain[i]] + 1];
-
- phi = 0;
- for(size_t j = 0; j < elementVector.size(); j++)
- {
- h = computeH(x, y, elementVector[j]);
- if(solution_type == 0)
- {
- g = computeG(x, y, elementVector[j], integrationType);
- }
- else
- {
- g = computeG_analytique(x, y, elementVector[j]);
- }
-
- phi = phi + h*elementVector[j].bcValue[0] - g*elementVector[j].bcValue[1];
- }
- vectorForNodeMap[i] = phi;
- }
- }
- for(std::size_t i = 0; i < NodeTags2d_domain.size(); i++)
- {
- nodeMapDomain[NodeTags2d_domain[i]] = vectorForNodeMap[i];
- }
-}
-
-void computeData(const std::vector<elementStruct> &elementVector, std::vector<size_t> &tags, std::vector<std::vector<double>> &data, const std::string integrationType, std::vector<int> &entities_tags_2d_domain_BEM, bool &solution_type)
-{
- double x = 0;
- double y = 0;
- double h = 0;
- double g = 0;
- double phi = 0;
- std::vector<double> barycenters;
- size_t elementCount = 0;
-
- for(size_t i = 0; i < entities_tags_2d_domain_BEM.size(); i++)
- {
- std::vector<int> elementTypes_for_entity;
- std::vector<std::vector<size_t> > elementTags_for_entity;
- std::vector<std::vector<size_t> > nodeTags_for_entity;
- gmsh::model::mesh::getElements(elementTypes_for_entity, elementTags_for_entity, nodeTags_for_entity, 2, entities_tags_2d_domain_BEM[i]);
-
- for(size_t l = 0; l < elementTypes_for_entity.size(); l++)
- {
- gmsh::model::mesh::getBarycenters(elementTypes_for_entity[l], entities_tags_2d_domain_BEM[i], false, false, barycenters);
-
-#pragma omp parallel for schedule(static) private(x,y,phi,h,g) shared(elementTags_for_entity,elementVector,tags,data,elementCount,barycenters,solution_type)
- for(size_t j = 0; j < elementTags_for_entity[l].size(); j++)
- {
- x = barycenters[3*j];
- y = barycenters[3*j+1];
-
- phi = 0;
- for(size_t k = 0; k < elementVector.size(); k++)
- {
- h = computeH(x, y, elementVector[k]);
- if(solution_type == 0)
- {
- g = computeG(x, y, elementVector[k], integrationType);
- }
- else
- {
- g = computeG_analytique(x, y, elementVector[k]);
- }
- phi += h*elementVector[k].bcValue[0] - g*elementVector[k].bcValue[1];
- }
- // Potential at the node is phi.
- tags[elementCount + j] = elementTags_for_entity[l][j];
- data[elementCount + j][0] = phi;
- }
- elementCount = elementCount + elementTags_for_entity[l].size();
- }
- }
-}
-
-void disp_element_node_data(std::vector<int> &nodes_commun_tag, std::vector<int> &elements_commun_tag)
-{
- std::cout << "problème dans la géométrie, noeuds communs = " << nodes_commun_tag.size() << std::endl;
- std::cout << "noeuds commun: "<< nodes_commun_tag.size() << ": ";
- for(size_t m = 0; m < nodes_commun_tag.size(); m++)
- {
- std::cout << nodes_commun_tag[m] << " ";
- }
- std::cout << std::endl;
- std::cout << "elements commun: " << elements_commun_tag.size() << ": ";
- for(size_t m = 0; m < elements_commun_tag.size(); m++)
- {
- std::cout << elements_commun_tag[m] << " ";
- }
- std::cout << std::endl;
- std::cout << std::endl;
-}
-
-void display_time(std::chrono::time_point<std::chrono::system_clock> &start, std::chrono::time_point<std::chrono::system_clock> &end, std::string &function_name, bool &show_time, const int domainTag)
-{
- //displays the time of the desired code section
- if(show_time)
- {
- std::cout << std::endl;
- std::cout << function_name << " in BEM domain " << domainTag << std::endl;
- std::cout << "Elapsed time in nanoseconds: "
- << std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count()
- << " ns" << std::endl;
-
- std::cout << "Elapsed time in microseconds: "
- << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
- << " µs" << std::endl;
-
- std::cout << "Elapsed time in milliseconds: "
- << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
- << " ms" << std::endl;
-
- std::cout << "Elapsed time in seconds: "
- << std::chrono::duration_cast<std::chrono::seconds>(end - start).count()
- << " sec";
- std::cout << std::endl;
- std::cout << std::endl;
- }
-}
-
-void convergenceFun(std::vector<int> &entities_tags_2d_domain_BEM, std::vector<size_t> &tags, std::vector<std::vector<double>> &phi)
-{
- double l = 1.0;
- double phiLeft = 100.0;
- double phiRight = 200.0;
-
- double error = 0.0;
-
- std::vector<double> barycenters;
-
- for(std::size_t i = 0; i < entities_tags_2d_domain_BEM.size(); i++)
- {
- std::vector<int> elementTypes_for_entity;
- std::vector<std::vector<size_t> > elementTags_for_entity;
- std::vector<std::vector<size_t> > nodeTags_for_entity;
- gmsh::model::mesh::getElements(elementTypes_for_entity, elementTags_for_entity, nodeTags_for_entity, 2, entities_tags_2d_domain_BEM[i]);
-
- for(std::size_t j = 0; j < elementTypes_for_entity.size(); j++)
- {
- gmsh::model::mesh::getBarycenters(elementTypes_for_entity[j], entities_tags_2d_domain_BEM[i], false, false, barycenters);
- for(std::size_t k = 0; k < elementTags_for_entity[i].size(); k++)
- {
- std::size_t tag = elementTags_for_entity[i][k];
- int index = 0;
- for(std::size_t l = 0; l < tags.size(); l++)
- {
- if(tag == tags[l])
- {
- index = l;
- break;
- }
- }
- double numericalPhi = phi[index][0];
-
- double x = barycenters[3*k];
- double exactPhi = phiLeft + (phiRight - phiLeft) * x / l;
-
- double relError = abs( (exactPhi - numericalPhi)/exactPhi );
- error += relError*relError;
- }
- }
- }
- error = sqrt(error/tags.size());
- std::cout << "\nNumber of elements = " << tags.size() << "\t-> Error = " << error << "\n\n";
-}
-
-double computeGradHx(const double &x, const double &y, const elementStruct &element, const std::string &integrationType)
-{
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
-
- double dx = xjp1-xj;
- double dy = yjp1-yj;
- double nx = dy;
- double ny = -dx;
- double norm = sqrt(nx*nx+ny*ny);
- nx = nx / norm;
- ny = ny / norm;
-
- // Computation of the integration points (parametric) and weights.
- std::vector<double> localCoord;
- std::vector<double> weights;
- gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
-
- // Computation of the sum (of integrand evaluated at integration point * weight)
- double gSum = 0;
- for(size_t k = 0; k < weights.size(); k++)
- {
- // Parametric coordinate
- double integration_variable = localCoord[3*k];
- // Global coordinates
- double hX = element.midX + element.deltaX*integration_variable/2.0;
- double hY = element.midY + element.deltaY*integration_variable/2.0;
-
- // Distance between mid point of element 'i' and integration point.
- double r = sqrt(pow(x - hX, 2.0) + pow(y - hY, 2.0));
- gSum += (-((1.0/2.0*pow(r, 2) - pow((xi - hX),2))/(pow(r, 4)))*nx + ((xi - hX)*(yi - hY)*ny)/(pow(r, 4))) * weights[k];
- }
-
- double gradHx = gSum * element.length / (2*M_PI);
-
- return gradHx;
-}
-double computeGradHy(const double &x, const double &y, const elementStruct &element, const std::string &integrationType)
-{
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
-
- double dx = xjp1-xj;
- double dy = yjp1-yj;
- double nx = dy;
- double ny = -dx;
- double norm = sqrt(nx*nx+ny*ny);
- nx = nx / norm;
- ny = ny / norm;
-
- // Computation of the integration points (parametric) and weights.
- std::vector<double> localCoord;
- std::vector<double> weights;
- gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
-
- // Computation of the sum (of integrand evaluated at integration point * weight)
- double gSum = 0;
- for(size_t k = 0; k < weights.size(); k++)
- {
- // Parametric coordinate
- double integration_variable = localCoord[3*k];
- // Global coordinates
- double hX = element.midX + element.deltaX*integration_variable/2.0;
- double hY = element.midY + element.deltaY*integration_variable/2.0;
-
- // Distance between mid point of element 'i' and integration point.
- double r = sqrt(pow(xi - hX, 2.0) + pow(yi - hY, 2.0));
- gSum += (-((1.0/2.0*pow(r, 2) - pow((yi - hY),2))/(pow(r, 4)))*ny + ((xi - hX)*(yi - hY)*nx)/(pow(r, 4))) * weights[k];
- }
-
- double gradHy = gSum * element.length / (2*M_PI);
-
- return gradHy;
-}
-double computeGradGx(const double &x, const double &y, const elementStruct &element, const std::string &integrationType)
-{
- // Computation of the integration points (parametric) and weights.
- std::vector<double> localCoord;
- std::vector<double> weights;
- gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
-
- // Computation of the sum (of integrand evaluated at integration point * weight)
- double gSum = 0;
- for(size_t k = 0; k < weights.size(); k++)
- {
- // Parametric coordinate
- double integration_variable = localCoord[3*k];
- // Global coordinates
- double gX = element.midX + element.deltaX*integration_variable/2.0;
- double gY = element.midY + element.deltaY*integration_variable/2.0;
-
- // Distance between mid point of element 'i' and integration point.
- double r = sqrt(pow(x - gX, 2.0) + pow(y - gY, 2.0));
- gSum += ((x-gX)/(pow(r, 2))) * weights[k];
- }
-
- double gradGx = gSum * element.length / (4*M_PI);
-
- return gradGx;
-}
-double computeGradGy(const double &x, const double &y, const elementStruct &element, const std::string &integrationType)
-{
- // Computation of the integration points (parametric) and weights.
- std::vector<double> localCoord;
- std::vector<double> weights;
- gmsh::model::mesh::getIntegrationPoints(1, integrationType, localCoord, weights);
-
- // Computation of the sum (of integrand evaluated at integration point * weight)
- double gSum = 0;
- for(size_t k = 0; k < weights.size(); k++)
- {
- // Parametric coordinate
- double integration_variable = localCoord[3*k];
- // Global coordinates
- double gX = element.midX + element.deltaX*integration_variable/2.0;
- double gY = element.midY + element.deltaY*integration_variable/2.0;
-
- // Distance between mid point of element 'i' and integration point.
- double r = sqrt(pow(x - gX, 2.0) + pow(y - gY, 2.0));
- gSum += ((y-gY)/(pow(r, 2))) * weights[k];
- }
-
- double gradGy = gSum * element.length / (4*M_PI);
-
- return gradGy;
-}
-
-double computeGradHx_analytique(const double &x, const double &y, const elementStruct &element)
-{
- double lj = element.length;
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
- double Pi = M_PI;
-
- double dx = xjp1-xj;
- double dy = yjp1-yj;
- double nx = dy;
- double ny = -dx;
- double norm = sqrt(nx*nx+ny*ny);
- nx = nx / norm;
- ny = ny / norm;
-
- double gradHx_analytique = (lj*(nx*(pow(xi,2) + xj*xjp1 - xi*(xj + xjp1) - (yi - yj)*(yi - yjp1)) -
- ny*(xj*yi + xjp1*yi - xjp1*yj - xj*yjp1 + xi*(-2*yi + yj + yjp1))))/
- (2.*Pi*(pow(xi - xj,2) + pow(yi - yj,2))*
- (pow(xi - xjp1,2) + pow(yi - yjp1,2)));
-
- return gradHx_analytique;
-}
-double computeGradHy_analytique(const double &x, const double &y, const elementStruct &element)
-{
- double lj = element.length;
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
- double Pi = M_PI;
-
- double dx = xjp1-xj;
- double dy = yjp1-yj;
- double nx = dy;
- double ny = -dx;
- double norm = sqrt(nx*nx+ny*ny);
- nx = nx / norm;
- ny = ny / norm;
-
- double gradHy_analytique = -(lj*(ny*(pow(xi,2) + xj*xjp1 - xi*(xj + xjp1) - (yi - yj)*(yi - yjp1)) +
- nx*(xj*yi + xjp1*yi - xjp1*yj - xj*yjp1 + xi*(-2*yi + yj + yjp1))))/
- (2.*Pi*(pow(xi - xj,2) + pow(yi - yj,2))*
- (pow(xi - xjp1,2) + pow(yi - yjp1,2)));
-
- return gradHy_analytique;
-}
-double computeGradGx_analytique(const double &x, const double &y, const elementStruct &element)
-{
- double lj = element.length;
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
- double Pi = M_PI;
-
- double gradGx_analytique = (lj*(-2*(yj - yjp1)*(atan(((xi - xj)*(xj - xjp1) + (yi - yj)*(yj - yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1))) +
- atan(((xi - xjp1)*(-xj + xjp1) + (yi - yjp1)*(-yj + yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1)))) -
- (xj - xjp1)*(log(pow(xi - xj,2) + pow(yi - yj,2)) - log(pow(xi - xjp1,2) + pow(yi - yjp1,2)))))/(4.*Pi*(pow(xj - xjp1,2) + pow(yj - yjp1,2)));
-
- return gradGx_analytique;
-}
-double computeGradGy_analytique(const double &x, const double &y, const elementStruct &element)
-{
- double lj = element.length;
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
- double Pi = M_PI;
-
- double gradGy_analytique = (lj*(2*(xj - xjp1)*(atan(((xi - xj)*(xj - xjp1) + (yi - yj)*(yj - yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1))) +
- atan(((xi - xjp1)*(-xj + xjp1) + (yi - yjp1)*(-yj + yjp1))/(xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1)))) -
- (yj - yjp1)*(log(pow(xi - xj,2) + pow(yi - yj,2)) - log(pow(xi - xjp1,2) + pow(yi - yjp1,2)))))/(4.*Pi*(pow(xj - xjp1,2) + pow(yj - yjp1,2)));
- return gradGy_analytique;
-}
-
-double computeG_analytique(const double &x, const double &y, const elementStruct &element)
-{
- double lj = element.length;
- double xi = x;
- double yi = y;
- double xj = element.x1;
- double xjp1 = element.x2;
- double yj = element.y1;
- double yjp1 = element.y2;
- double Pi = M_PI;
-
- double g_analytique = (lj*(-2*(pow(xj - xjp1,2) + pow(yj - yjp1,2)) + 2*(-(xj*yi) + xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*yjp1)* atan(((xi - xjp1)*(xj - xjp1) + (yi - yjp1)*(yj - yjp1))/(-(xj*yi) + xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*yjp1)) + 2*(-(xjp1*yi) - xi*yj + xjp1*yj + xj*(yi - yjp1) + xi*yjp1)* atan(((xi - xj)*(xj - xjp1) + (yi - yj)*(yj - yjp1))/ (xjp1*yi + xi*yj - xjp1*yj - xi*yjp1 + xj*(-yi + yjp1))) + (-((xi - xj)*(xj - xjp1)) - (yi - yj)*(yj - yjp1))* log(pow(xi - xj,2) + pow(yi - yj,2)) + ((xi - xjp1)*(xj - xjp1) + (yi - yjp1)*(yj - yjp1))* log(pow(xi - xjp1,2) + pow(yi - yjp1,2))))/(4.*Pi*(pow(xj - xjp1,2) + pow(yj - yjp1,2)));
-
- return g_analytique;
-}
-
-void electricField(int &nb2dElements, std::vector<int> &entities_tags_2d_domain_BEM, const std::vector<elementStruct> &elementVector, std::vector<size_t> &tags, const std::string &integrationType, std::vector<std::vector<double>> &data_Electric_field, bool &solution_type)
-{
- double x = 0;
- double y = 0;
- double grad_phi_x = 0;
- double grad_phi_y = 0;
- double h_x = 0;
- double g_x = 0;
- double h_y = 0;
- double g_y = 0;
- std::vector<double> barycenters;
- size_t elementCount = 0;
-
- for(size_t i = 0; i < entities_tags_2d_domain_BEM.size(); i++)
- {
- std::vector<int> elementTypes_for_entity;
- std::vector<std::vector<size_t> > elementTags_for_entity;
- std::vector<std::vector<size_t> > nodeTags_for_entity;
- gmsh::model::mesh::getElements(elementTypes_for_entity, elementTags_for_entity, nodeTags_for_entity, 2, entities_tags_2d_domain_BEM[i]);
-
- for(size_t l = 0; l < elementTypes_for_entity.size(); l++)
- {
- gmsh::model::mesh::getBarycenters(elementTypes_for_entity[l], entities_tags_2d_domain_BEM[i], false, false, barycenters);
-#pragma omp parallel for schedule(static) private(x,y,grad_phi_x,grad_phi_y,h_x,g_x,h_y,g_y) shared(elementTags_for_entity,elementVector,tags,data_Electric_field,elementCount,barycenters,solution_type,integrationType)
- for(size_t j = 0; j < elementTags_for_entity[l].size(); j++)
- {
- x = barycenters[3*j];
- y = barycenters[3*j+1];
-
- grad_phi_x = 0;
- grad_phi_y = 0;
- for(size_t k = 0; k < elementVector.size(); k++)
- {
- if(solution_type == 0)
- {
- h_x = computeGradHx(x, y, elementVector[k], integrationType);
- g_x = computeGradGx(x, y, elementVector[k], integrationType);
- h_y = computeGradHy(x, y, elementVector[k], integrationType);
- g_y = computeGradGy(x, y, elementVector[k], integrationType);
- }
- else
- {
- h_x = computeGradHx_analytique(x, y, elementVector[k]);
- g_x = computeGradGx_analytique(x, y, elementVector[k]);
- h_y = computeGradHy_analytique(x, y, elementVector[k]);
- g_y = computeGradGy_analytique(x, y, elementVector[k]);
- }
- grad_phi_x = grad_phi_x + h_x*elementVector[k].bcValue[0] - g_x*elementVector[k].bcValue[1];
- grad_phi_y = grad_phi_y + h_y*elementVector[k].bcValue[0] - g_y*elementVector[k].bcValue[1];
- }
- tags[elementCount + j] = elementTags_for_entity[l][j];
- data_Electric_field[elementCount + j][0] = -grad_phi_x;
- data_Electric_field[elementCount + j][1] = -grad_phi_y;
- data_Electric_field[elementCount + j][2] = 0;
- }
- elementCount = elementCount + elementTags_for_entity[l].size();
- }
- }
-}
diff --git a/srcs/BEM/functionsBEM.hpp b/srcs/BEM/functionsBEM.hpp
index ea63482973bfa134d5c6af3b28854631f1d85b2a..c2cce535a33e45229b983a63ae0000941a7ff197 100644
--- a/srcs/BEM/functionsBEM.hpp
+++ b/srcs/BEM/functionsBEM.hpp
@@ -11,41 +11,199 @@
// Each element has all its information gathered as in this structure.
struct elementStruct{
- int tag; // tag of the element
- std::string bcType; // "dirichlet" or "neumann" (at the end, should replace this by an int or boolean)
+ std::size_t tag; // tag of the element
double bcValue [2]; // values of [u, u_n]
size_t nodeTags [2]; // tags of the [first, second] nodes of the element
- size_t midNodeTags;
+ size_t midNodeTag;
double length;
double x1, y1, x2, y2;
double midX, midY;
double deltaX, deltaY;
};
+/*
+ * Main
+ */
-void solverBEM(std::map<int, double> &electrostaticPressure, int & nbViews, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, bool postProcessing, const int iteration, bool untangleMesh);
-void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int & nbViews, bool postProcessing, const int iteration, const int domainTag);
-void dispNodeCoordFun(const std::vector<size_t> &allNodeTags, const std::vector<double> &allCoord);
-void dispMatricesFun(const Eigen::MatrixXd &A, const Eigen::MatrixXd &B, const Eigen::MatrixXd &c, const Eigen::MatrixXd &b, const Eigen::MatrixXd &x);
-void dispValuesFun(const std::vector<elementStruct> &elementVector);
-void getBC(std::map<size_t, int> &allNodeIndex, const std::vector<double> &allCoord, std::vector<elementStruct> &elementVector, int &dirichletCount, const bool dispHierarchy, std::vector<size_t> &NodeTags2d_domain, std::vector<double> &NodeCoord2d_domain, std::vector<int> &elementTypes2d_domain_BEM, std::vector<std::vector<size_t>> &elementTags2d_domain_BEM, std::vector<std::vector<size_t>> &elnodeTags2d_domain_BEM, int &entities_tags_for_physical_group_domain_BEM, std::vector<int> &entities_tags_2d_domain_BEM, std::map<int, std::pair<double, double>> &boundaryDisplacementMap, std::map<std::size_t, std::pair<bool,bool>> &nodeOnBoundary, const int domainTag);
-void computeGeometricParam(elementStruct &element, std::map<std::size_t, int> &allNodeIndex, const std::vector<double> &allCoord, std::pair<double, double> &nodalDisp1, std::pair<double, double> &nodalDisp2);
-double computeH(const double &x, const double &y, const elementStruct &element);
-double computeG(const double &x, const double &y, const elementStruct &element, const std::string &integrationType);
-void computeData(const std::vector<elementStruct> &elementVector, std::vector<size_t> &tags, std::vector<std::vector<double>> &data, const std::string integrationType, std::vector<int> &entities_tags_2d_domain_BEM, bool &solution_type);
-void displayData(const std::string viewName, const std::string modelName, const std::string fileName, const std::vector<size_t> &tags, const std::vector<std::vector<double>> &data);
-void compute_internal_phi(std::vector<double> &allCoord, std::vector<size_t> &NodeTags2d_domain, std::map<size_t, int> &allNodeIndex, std::vector<elementStruct> &elementVector, std::map<int,double> &nodeMapDomain, const std::string integrationType, std::map<std::size_t, std::pair<bool,bool>> &nodeOnBoundary, bool &solution_type);
-void element_node_data(std::vector<std::vector<std::vector<double>>> &data, std::vector<std::vector<size_t>> &elementTags2d_domain_BEM, std::vector<int> &elementTypes2d_domain_BEM, std::vector<std::vector<size_t>> &nodeTags2d_domain_BEM, std::vector<elementStruct> &elementVector, std::map<int, double> &nodeMapDomain, bool &show_element_node_data, int &entities_tags_for_physical_group_domain_BEM, std::map<std::size_t, std::pair<bool,bool>> &nodeOnBoundary);
-void disp_element_node_data(std::vector<int> &nodes_commun_tag, std::vector<int> &elements_commun_tag);
-void display_time(std::chrono::time_point<std::chrono::system_clock> &start, std::chrono::time_point<std::chrono::system_clock> &end, std::string &function_name, bool &show_time, const int domainTag);
-void convergenceFun(std::vector<int> &entities_tags_2d_domain_BEM, std::vector<size_t> &tags, std::vector<std::vector<double>> &phi);
-double computeGradHx(const double &x, const double &y, const elementStruct &element, const std::string &integrationType);
-double computeGradHy(const double &x, const double &y, const elementStruct &element, const std::string &integrationType);
-double computeGradGx(const double &x, const double &y, const elementStruct &element, const std::string &integrationType);
-double computeGradGy(const double &x, const double &y, const elementStruct &element, const std::string &integrationType);
-double computeGradHx_analytique(const double &x, const double &y, const elementStruct &element);
-double computeGradHy_analytique(const double &x, const double &y, const elementStruct &element);
-double computeGradGx_analytique(const double &x, const double &y, const elementStruct &element);
-double computeGradGy_analytique(const double &x, const double &y, const elementStruct &element);
-double computeG_analytique(const double &x, const double &y, const elementStruct &element);
-void electricField(int &nb2dElements, std::vector<int> &entities_tags_2d_domain_BEM, const std::vector<elementStruct> &elementVector, std::vector<size_t> &tags, const std::string &integrationType, std::vector<std::vector<double>> &data_Electric_field, bool &solution_type);
+// Fills the "electrostaticPressure" map that gives the electrostatic pressure associated with each (boundary) element tag.
+// "nbViews" is the counter of views (useful if other views have to be generated than with this solver). "nodalDisplacements" is
+// a map associating the (x,y) displacement to a serie of node tags (if the position of the nodes are not striclty the one
+// extracted from the .geo file). "postProcessing" indicates whether the post-processing should be done ('true') or not
+// ('false'). "iteration" represents the number of times the solver has been called (some information is printed on the first
+// iteration). "untangleMesh" unstangles every BEM domain if set to 'true' (in case the "nodalDisplacements" map contains
+// non-zero displacements).
+void solverBEM(std::map<int, double> &electrostaticPressure,
+ int & nbViews,
+ std::map<int,std::pair<double,double>> & nodalDisplacements,
+ const bool postProcessing,
+ const int iteration,
+ const bool untangleMesh);
+
+// Fills the "electrostaticPressure" map. The "nodalDisplacements", "postProcessing" and "iteration" arguments are described in
+// the description of the function "solverBEM". "domainName" is used for extracting everything in the appropriate independant BEM
+// domain, "phiViewTag" and "electricFieldViewTag" are the view tag of the potential and of the electric field respectively.
+// "computePhi" and "computeElectricField" indicates whether the potential and/or the electric field should be computed ('true')
+// or not ('false').
+void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
+ const int nbViews,
+ std::map<int,std::pair<double,double>> & nodalDisplacements,
+ const bool postProcessing,
+ const int iteration,
+ const std::string domainName,
+ const int phiViewTag,
+ const int electricFieldViewTag,
+ const bool computePhi,
+ const bool computeElectricField);
+
+/*
+ * Boundary element
+ */
+
+// Fills the "boundaryConditions" map associating to each physical group (a string) the type of boundary condition ('true':
+// dirichlet, 'false': neumann) and its value (if any). "domainName" is the name of the current BEM domain.
+void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
+ const std::string domainName);
+
+// Fills the "elementVector" vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the
+// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also fills the
+// "physicalGroups" map associating to each (physical group) name, its tag and dimension. "domainName" is the name of the current
+// BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement to a serie of node tags. Returns the number of
+// 'dirichlet elements' (useful for the construction/resolution of the system).
+std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
+ std::map<std::string, std::pair<int, int>> &physicalGroups,
+ std::map<int, std::pair<double, double>> &nodalDisplacements,
+ const std::string domainName);
+
+// Computes the useful x of a given element and stores it in "element". "nodeIndinces" is a map associating the index of any
+// node tag for accessing its coordinates in the "nodeCoords" vector. "nodalDisp1" and "nodalDisp2" are pairs representing the
+// potential displacement of the first and second node (respectivelly) of the element. If an element has 3 nodes, the first and
+// second nodes are considered to be the extremities of that element.
+void computeGeometricParam(elementStruct &element,
+ std::map<std::size_t, std::size_t> &nodeIndices,
+ const std::vector<double> &nodeCoords,
+ const std::pair<double, double> &nodalDisp1,
+ const std::pair<double, double> &nodalDisp2);
+
+// Fills the "A", "B", "c" matrices (/vector) that represent the linear system. "elementVector" is a vector storing the boundary
+// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "nbDirichletElements" is the number of
+// 'dirichlet elements'. "solutionType" indicates whether the elements of the matrices are computed analytically ('true') or
+// numerically ('false'). "integrationType" is the type of numerical integration used (in case of numerical computation).
+void fillSystem(Eigen::MatrixXd &A,
+ Eigen::MatrixXd &B,
+ Eigen::MatrixXd &c,
+ const std::vector<elementStruct> &elementVector,
+ const std::size_t nbDirichletElements,
+ const std::string integrationType,
+ const bool solutionType);
+
+// Fills the "electrostaticPressure" map associating each (boundary) element tag with its electrostatic pressure. Also sets the
+// unknown of each element (either the potential of the electric field) of the "elementVector" vector storing the boundary
+// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "x" is the solution of the linear
+// system. "nbDirichletElements" is the number of 'dirichlet elements', "domainName" is the name of the considered BEM domain.
+void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
+ std::vector<elementStruct> &elementVector,
+ const Eigen::MatrixXd &x,
+ const std::size_t nbDirichletElements,
+ const std::string domainName);
+
+/*
+ * Post processing
+ */
+
+int getModel(std::vector<size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ std::vector<double> &domainNodeCoords,
+ std::vector<int> &domainEntityTags,
+ std::map<int, std::size_t> &domainEntityIndices,
+ std::vector<std::vector<int>> &domainElementTypes,
+ std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ std::map<std::string, std::pair<int, int>> &physicalGroups,
+ const std::string domainName);
+
+void computeInternalPhi(std::map<int, double> &internalNodalPhi,
+ const std::vector<double> &domainNodeCoords,
+ const std::vector<size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ const std::vector<elementStruct> &elementVector,
+ const std::string integrationType,
+ std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
+ const bool &solutionType);
+
+void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
+ const std::vector<elementStruct> &elementVector);
+
+void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ const std::vector<elementStruct> &elementVector,
+ const int domainPhysicalGroupTag);
+
+void phiInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ std::map<std::size_t, std::vector<std::size_t>> &position,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags);
+
+void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ std::map<int, double> &internalNodalPhi);
+
+void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ const std::vector<elementStruct> &elementVector,
+ std::map<int, double> &internalNodalPhi,
+ const int domainPhysicalGroupTag,
+ std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements);
+
+void computeElementData(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ std::vector<std::vector<std::vector<std::vector<double>>>> &electricField,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<int> &domainEntityTags,
+ const std::vector<elementStruct> &elementVector,
+ const bool solutionType,
+ const std::string integrationType,
+ const bool computePhi,
+ const bool computeElectricField);
+
+void displayData(const int viewTag,
+ const std::string modelName,
+ const std::string fileName,
+ const std::string dataType,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &data,
+ const int viewIndex);
+
+/*
+ * Verification
+ */
+
+void displayTime(const double start,
+ const double end,
+ std::string functionName, const bool showTime);
+
+void convergenceFun(const std::vector<int> &domainEntityTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalPhi,
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalElectricField);
+
+
+/*
+ * Computations
+ */
+
+double computeAnalyticalH(const double x, const double y, const elementStruct &element);
+double computeNumericalG(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradHx(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradHy(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradGx(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradGy(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeAnalyticalGradHx(const double x, const double y, const elementStruct &element);
+double computeAnalyticalGradHy(const double x, const double y, const elementStruct &element);
+double computeAnalyticalGradGx(const double x, const double y, const elementStruct &element);
+double computeAnalyticalGradGy(const double x, const double y, const elementStruct &element);
+double computeAnalyticalG(const double x, const double y, const elementStruct &element);
diff --git a/srcs/BEM/mainBEM.cpp b/srcs/BEM/mainBEM.cpp
index 9f35b0d70edee81bf0e49e41cb5b9ba0c078fbe1..02035781a5f8ac07e30bc7e0ac8848a8994c4116 100644
--- a/srcs/BEM/mainBEM.cpp
+++ b/srcs/BEM/mainBEM.cpp
@@ -18,25 +18,15 @@ int main(int argc, char **argv)
// will map the 1d element tag onto the corresponding electroPressure resulting from BEM code
std::map<int, double> electrostaticPressure;
-
int nbViews = 0;
+ std::map<int,std::pair<double,double>> nodalDisplacements;
+ const bool postProcessing = true;
+ const bool meshUntangler = false;
- std::map<int,std::pair<double,double>> boundaryDisplacementMap;
-
- //solverBEM(argc, argv, finalElementTags, electrostaticPressure, nbViews);
- solverBEM(electrostaticPressure, nbViews, boundaryDisplacementMap, true, 1, false); // iteration number randomly set to 1
-
- std::map<int, double>::iterator it;
-
- /*for (it = electrostaticPressure.begin(); it != electrostaticPressure.end(); it++)
- {
- std::cout << it->first // string (key)
- << ':'
- << it->second // string's value
- << std::endl;
- }*/
+ solverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, 1, meshUntangler); // iteration number randomly set to 1
- gmsh::fltk::run();
+ if(postProcessing)
+ gmsh::fltk::run();
gmsh::finalize();
return 0;
diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c5854559391eff43301a958a05c5861d9c752c03
--- /dev/null
+++ b/srcs/BEM/postProcessing.cpp
@@ -0,0 +1,468 @@
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+#include <sstream>
+#include <Eigen/Dense>
+#include <cmath>
+#include <iomanip>
+#include <stdio.h>
+#include <chrono>
+#include "functionsBEM.hpp"
+#include <omp.h>
+#include <algorithm>
+#include <list>
+
+int getModel(std::vector<std::size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ std::vector<double> &domainNodeCoords,
+ std::vector<int> &domainEntityTags,
+ std::map<int, std::size_t> &domainEntityIndices,
+ std::vector<std::vector<int>> &domainElementTypes,
+ std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ std::map<std::string, std::pair<int, int>> &physicalGroups,
+ const std::string domainName)
+{
+ // Retrieves everything only for the domain of interrest.
+ int physicalGroupDim = physicalGroups[domainName].first;
+ int physicalGroupTag = physicalGroups[domainName].second;
+
+ gmsh::model::mesh::getNodesForPhysicalGroup(physicalGroupDim, physicalGroupTag, domainNodeTags, domainNodeCoords);
+
+ // The "domainNodeIndices" map is used for finding the position of the coordinates of a node in "domainNodeCoords"
+ // from its tag.
+ for(std::size_t index = 0; index < domainNodeTags.size(); index++)
+ domainNodeIndices[domainNodeTags[index]] = index;
+
+ gmsh::model::getEntitiesForPhysicalGroup(physicalGroupDim, physicalGroupTag, domainEntityTags);
+ domainElementTypes.resize(domainEntityTags.size());
+ domainElementTags.resize(domainEntityTags.size());
+ for(std::size_t j = 0; j < domainEntityTags.size(); j++)
+ {
+ domainEntityIndices[domainEntityTags[j]] = j;
+
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> nodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, physicalGroupDim, domainEntityTags[j]);
+
+ domainElementTypes[j].resize(elementTypes.size());
+ domainElementTags[j].resize(elementTypes.size());
+ for(std::size_t k = 0; k < elementTypes.size(); k++)
+ {
+ int nbNodes = nodeTags[k].size() / elementTags[k].size();
+ std::vector<std::size_t> elementNodeTags(nbNodes);
+
+ domainElementTypes[j][k] = elementTypes[k];
+ domainElementTags[j][k].resize(elementTags[k].size());
+
+ for(std::size_t l = 0; l < elementTags[k].size(); l++)
+ {
+ domainElementTags[j][k][l] = elementTags[k][l];
+ int n = 0;
+ for(std::size_t m = l*nbNodes; m < (l+1)*nbNodes; m++)
+ {
+ elementNodeTags[n] = nodeTags[k][m];
+ n++;
+ }
+
+ domainElementNodeTags[elementTags[k][l]] = elementNodeTags;
+ }
+ }
+ }
+
+ return physicalGroupTag;
+}
+
+
+void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
+ const std::vector<elementStruct> &elementVector)
+{
+ for(std::size_t i = 0; i < elementVector.size(); i++)
+ {
+ if(!boundaryNodeToElements.count(elementVector[i].nodeTags[0]))
+ boundaryNodeToElements[elementVector[i].nodeTags[0]] = {i};
+ else
+ boundaryNodeToElements[elementVector[i].nodeTags[0]].push_back(i);
+
+ if(!boundaryNodeToElements.count(elementVector[i].nodeTags[1]))
+ boundaryNodeToElements[elementVector[i].nodeTags[1]] = {i};
+ else
+ boundaryNodeToElements[elementVector[i].nodeTags[1]].push_back(i);
+
+ if(elementVector[i].midNodeTag != 0)
+ boundaryNodeToElements[elementVector[i].midNodeTag] = {i};
+ }
+}
+
+void computeInternalPhi(std::map<int, double> &internalNodalPhi,
+ const std::vector<double> &domainNodeCoords,
+ const std::vector<std::size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ const std::vector<elementStruct> &elementVector,
+ const std::string integrationType,
+ std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
+ const bool &solutionType)
+{
+ double x = 0;
+ double y = 0;
+ double phi = 0;
+ double h = 0;
+ double g = 0;
+ std::vector<double> nodeVectorDomain(domainNodeTags.size(), 0);
+
+ #pragma omp parallel for schedule(static) private(x, y, phi, h, g) shared(domainNodeTags, domainNodeCoords, elementVector, internalNodalPhi, boundaryNodeToElements, solutionType)
+ for(std::size_t i = 0; i < domainNodeTags.size(); i++)
+ {
+ if(!boundaryNodeToElements.count(domainNodeTags[i]))
+ {
+ // Coordinates of the node
+ x = domainNodeCoords[3*domainNodeIndices[domainNodeTags[i]]];
+ y = domainNodeCoords[3*domainNodeIndices[domainNodeTags[i]] + 1];
+
+ phi = 0;
+ for(size_t j = 0; j < elementVector.size(); j++)
+ {
+ h = computeAnalyticalH(x, y, elementVector[j]);
+ if(!solutionType)
+ g = computeNumericalG(x, y, elementVector[j], integrationType);
+ else
+ g = computeAnalyticalG(x, y, elementVector[j]);
+
+ phi += h*elementVector[j].bcValue[0] - g*elementVector[j].bcValue[1];
+ }
+ nodeVectorDomain[i] = phi;
+ }
+ }
+
+ for(std::size_t i = 0; i < domainNodeTags.size(); i++)
+ internalNodalPhi[domainNodeTags[i]] = nodeVectorDomain[i];
+}
+
+void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ const std::vector<elementStruct> &elementVector,
+ std::map<int, double> &internalNodalPhi,
+ const int domainPhysicalGroupTag,
+ std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements)
+{
+ // This map may be replaced in the future. It allows to retain the position of the elements of the elementTag2d loop and linked their position to the loop of the elementTag2d_boundary (the tags of the elements in common with the boundary).
+ // This map is the link between the array of tags of the 2d elements of the boundary and that of the tags of all the 2d elements.
+ std::map<std::size_t, std::vector<std::size_t>> position;
+ phiInitialization(phi, position, domainElementTypes, domainElementTags);
+
+ // // The first part of the algorithm consists in recovering the tag of the 2d elements of the domain which have at least one node in common with the boundary
+ // // We initialise the first dimension of the vector that will contain all 2d elements in common with the boundary with the number of element types present in the domain.
+ std::vector<size_t> domainBoundaryElementTags;
+ boundaryElementTags(domainBoundaryElementTags, domainElementTypes, domainElementTags, elementVector, domainPhysicalGroupTag);
+
+ // The second part of the algorithm will consist in assigning the correct potential value to the nodes of the elements in common with the boundary.
+ // The first loop is on the elementType in 2d in common with the boundary.
+ #pragma omp parallel for schedule(guided)
+ for(std::size_t i = 0; i < domainBoundaryElementTags.size(); i++)
+ {
+ std::size_t elementTag = domainBoundaryElementTags[i];
+ int entityTagIndex = position[elementTag][0];
+ int elementTypeIndex = position[elementTag][1];
+ int elementTagIndex = position[elementTag][2];
+ int nbNodes = phi[entityTagIndex][elementTypeIndex][elementTagIndex].size();
+ std::vector<std::size_t> nodeTags = domainElementNodeTags[elementTag];
+
+ // The third loop is on the nodes of a specific elementTag in 2d in common with the boundary.
+ for(int j = 0; j < nbNodes; j++)
+ {
+ std::size_t elementNodeTag = nodeTags[j];
+
+ if(boundaryNodeToElements.count(elementNodeTag))
+ {
+ std::vector<std::size_t> commonBoundaryElementIndices = boundaryNodeToElements[elementNodeTag];
+
+ // Mid node of a segment on the boundary => take the value of the element.
+ if(commonBoundaryElementIndices.size() == 1)
+ {
+ phi[entityTagIndex][elementTypeIndex][elementTagIndex][j] = elementVector[commonBoundaryElementIndices[0]].bcValue[0];
+ }
+ else if(commonBoundaryElementIndices.size() == 2)
+ {
+ double currentPhi = 0.0;
+ int count = 0;
+
+ // First element.
+ int adjacentElementIndex = commonBoundaryElementIndices[0];
+ std::size_t adjacentNodeTag = elementVector[adjacentElementIndex].nodeTags[0];
+ if(adjacentNodeTag == elementNodeTag)
+ {
+ adjacentNodeTag = elementVector[adjacentElementIndex].nodeTags[1];
+ }
+
+ // if 'adjacentNodeTag' is in 'nodeTags'
+ if(std::count(nodeTags.begin(), nodeTags.end(), adjacentNodeTag))
+ {
+ count++;
+ }
+ currentPhi += elementVector[adjacentElementIndex].bcValue[0];
+
+
+ // Second element.
+ adjacentElementIndex = commonBoundaryElementIndices[1];
+ adjacentNodeTag = elementVector[adjacentElementIndex].nodeTags[0];
+ if(adjacentNodeTag == elementNodeTag)
+ {
+ adjacentNodeTag = elementVector[adjacentElementIndex].nodeTags[1];
+ }
+
+ // if 'adjacentNodeTag' is in 'nodeTags'
+ if(std::count(nodeTags.begin(), nodeTags.end(), adjacentNodeTag))
+ {
+ if(count == 0)
+ {
+ currentPhi = elementVector[adjacentElementIndex].bcValue[0];
+ }
+ else
+ {
+ currentPhi += elementVector[adjacentElementIndex].bcValue[0];
+ currentPhi /= 2;
+ }
+ }
+ else
+ {
+ if(count == 0)
+ {
+ currentPhi += elementVector[adjacentElementIndex].bcValue[0];
+ currentPhi /= 2;
+ }
+ }
+
+ phi[entityTagIndex][elementTypeIndex][elementTagIndex][j] = currentPhi;
+ }
+ }
+ }
+ }
+
+ //The last step is to loop over all the 2d elements in the domain. If the HUGE_VAL value is still present, this means that the element node does not belong to the boundary (since the value of the potential has not been modified) and is therefore part of the domain. We then use the map which associates the value of the potential with the tag of the node inside the domain.
+ boundaryElementFinalization(phi, domainElementTags, domainElementNodeTags, internalNodalPhi);
+}
+
+
+void phiInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ std::map<std::size_t, std::vector<std::size_t>> &position,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags)
+{
+ phi.resize(domainElementTypes.size());
+ for(std::size_t i = 0; i < domainElementTypes.size(); i++)
+ {
+ phi[i].resize(domainElementTypes[i].size());
+
+ #pragma omp parallel for schedule(guided) shared(phi, domainElementTypes, domainElementTags, position)
+ for(std::size_t j = 0; j < domainElementTypes[i].size(); j++)
+ {
+ std::string elementName;
+ int dim;
+ int order;
+ int nbNodes;
+ std::vector<double> localNodeCoord;
+ int nbPrimaryNodes;
+ gmsh::model::mesh::getElementProperties(domainElementTypes[i][j], elementName, dim, order, nbNodes, localNodeCoord, nbPrimaryNodes);
+
+ phi[i][j].resize(domainElementTags[i][j].size());
+ for(std::size_t k = 0; k < domainElementTags[i][j].size(); k++)
+ {
+ phi[i][j][k].resize(nbNodes, HUGE_VAL);
+ std::vector<std::size_t> pos = {i, j, k};
+
+ #pragma omp critical
+ position[domainElementTags[i][j][k]] = pos;
+ }
+ }
+ }
+}
+
+void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ const std::vector<elementStruct> &elementVector,
+ const int domainPhysicalGroupTag)
+{
+ #pragma omp parallel
+ {
+ std::vector<std::size_t> localElementTags;
+
+ // This loop goes through all the nodes of the boundary (node 1 of each element). We then use getElementsByCoordinates to obtain all the elements in common with the node considered. We only take the 2d elements.
+ #pragma omp for schedule(guided)
+ for(std::size_t i = 0; i < elementVector.size(); i++)
+ {
+ std::vector<size_t> elementTags;
+ gmsh::model::mesh::getElementsByCoordinates(elementVector[i].x1, elementVector[i].y1, 0, elementTags, 2);
+
+ // We loop on the 2d elements which are in common with the considered node of the boundary. Depending on their type, they are added to the vector comprising all the 2d elements in contact with the boundary.
+ for(std::size_t j = 0; j < elementTags.size(); j++)
+ {
+ int elementType;
+ std::vector<size_t> nodeTags;
+ int entityDim;
+ int entityTag;
+ std::vector<int> physicalTags;
+ #pragma omp critical
+ {
+ gmsh::model::mesh::getElement(elementTags[j], elementType, nodeTags, entityDim, entityTag);
+ gmsh::model::getPhysicalGroupsForEntity(entityDim, entityTag, physicalTags);
+ }
+
+ for(std::size_t k = 0; k < physicalTags.size(); k++)
+ if(physicalTags[k] == domainPhysicalGroupTag)
+ localElementTags.push_back(elementTags[j]);
+ }
+ }
+
+ #pragma omp critical
+ domainBoundaryElementTags.insert(domainBoundaryElementTags.end(), localElementTags.begin(), localElementTags.end());
+ }
+
+
+ // Allows duplicate tags to be removed. If an element has several nodes with the border, it has been counted several times before.
+ sort( domainBoundaryElementTags.begin(), domainBoundaryElementTags.end() );
+ domainBoundaryElementTags.erase( unique( domainBoundaryElementTags.begin(), domainBoundaryElementTags.end() ), domainBoundaryElementTags.end() );
+}
+
+
+void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ std::map<int, double> &internalNodalPhi)
+{
+ for(std::size_t i = 0; i < domainElementTags.size(); i++)
+ {
+ for(std::size_t j = 0; j < domainElementTags[i].size(); j++)
+ {
+ #pragma omp parallel for schedule(guided)
+ for(std::size_t k = 0; k < domainElementTags[i][j].size(); k++)
+ {
+ std::vector<std::size_t> nodeTags = domainElementNodeTags[domainElementTags[i][j][k]];
+ for(std::size_t l = 0; l < nodeTags.size(); l++)
+ {
+ std::size_t elementNodeTag = nodeTags[l];
+ if(phi[i][j][k][l] == HUGE_VAL && internalNodalPhi.count(elementNodeTag))
+ phi[i][j][k][l] = internalNodalPhi[elementNodeTag];
+ }
+ }
+ }
+ }
+}
+
+
+void computeElementData(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ std::vector<std::vector<std::vector<std::vector<double>>>> &electricField,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<int> &domainEntityTags,
+ const std::vector<elementStruct> &elementVector,
+ const bool solutionType,
+ const std::string integrationType,
+ const bool computePhi,
+ const bool computeElectricField)
+{
+ std::vector<double> barycenters;
+
+ phi.resize(domainEntityTags.size());
+ electricField.resize(domainEntityTags.size());
+ for(std::size_t i = 0; i < domainEntityTags.size(); i++)
+ {
+ phi[i].resize(domainElementTypes[i].size());
+ electricField[i].resize(domainElementTypes[i].size());
+ for(std::size_t j = 0; j < domainElementTypes[i].size(); j++)
+ {
+ phi[i][j].resize(domainElementTags[i][j].size());
+ electricField[i][j].resize(domainElementTags[i][j].size());
+ gmsh::model::mesh::getBarycenters(domainElementTypes[i][j], domainEntityTags[i], false, false, barycenters);
+ for(size_t k = 0; k < domainElementTags[i][j].size(); k++)
+ {
+ double x = barycenters[3*k];
+ double y = barycenters[3*k + 1];
+
+ double currentPhi = 0;
+ double gradPhiX = 0;
+ double gradPhiY = 0;
+
+ #pragma omp parallel shared(currentPhi, gradPhiX, gradPhiY)
+ {
+ double localCurrentPhi = 0;
+ double localGradPhiX = 0;
+ double localGradPhiY = 0;
+
+ double h, hX, hY;
+ double g, gX, gY;
+
+ #pragma omp for schedule(guided)
+ for(std::size_t el = 0; el < elementVector.size(); el++)
+ {
+ if(computePhi)
+ {
+ h = computeAnalyticalH(x, y, elementVector[el]);
+ if(!solutionType)
+ g = computeNumericalG(x, y, elementVector[el], integrationType);
+ else
+ g = computeAnalyticalG(x, y, elementVector[el]);
+
+ localCurrentPhi += h*elementVector[el].bcValue[0] - g*elementVector[el].bcValue[1];
+ }
+ if(computeElectricField)
+ {
+ if(!solutionType)
+ {
+ hX = computeNumericalGradHx(x, y, elementVector[el], integrationType);
+ gX = computeNumericalGradGx(x, y, elementVector[el], integrationType);
+ hY = computeNumericalGradHy(x, y, elementVector[el], integrationType);
+ gY = computeNumericalGradGy(x, y, elementVector[el], integrationType);
+ }
+ else
+ {
+ hX = computeAnalyticalGradHx(x, y, elementVector[el]);
+ gX = computeAnalyticalGradGx(x, y, elementVector[el]);
+ hY = computeAnalyticalGradHy(x, y, elementVector[el]);
+ gY = computeAnalyticalGradGy(x, y, elementVector[el]);
+ }
+
+ localGradPhiX += hX*elementVector[el].bcValue[0] - gX*elementVector[el].bcValue[1];
+ localGradPhiY += hY*elementVector[el].bcValue[0] - gY*elementVector[el].bcValue[1];
+ }
+ }
+
+ #pragma omp critical
+ {
+ currentPhi += localCurrentPhi;
+ gradPhiX += localGradPhiX;
+ gradPhiY += localGradPhiY;
+ }
+ }
+
+ phi[i][j][k].resize(1, currentPhi);
+ electricField[i][j][k].resize(3);
+ electricField[i][j][k][0] = -gradPhiX;
+ electricField[i][j][k][1] = -gradPhiY;
+ electricField[i][j][k][2] = 0.0;
+ }
+ }
+ }
+}
+
+void displayData(const int viewTag,
+ const std::string modelName,
+ const std::string fileName,
+ const std::string dataType,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &data,
+ const int viewIndex)
+{
+ for(std::size_t i = 0; i < data.size(); i++)
+ for(std::size_t j = 0; j < data[i].size(); j++)
+ gmsh::view::addModelData(viewTag, 0 , modelName, dataType, domainElementTags[i][j], data[i][j]);
+
+ gmsh::view::write(viewTag, fileName);
+
+ gmsh::option::setNumber("View[" + std::to_string(viewIndex) + "].IntervalsType", 3);
+ gmsh::option::setNumber("View[" + std::to_string(viewIndex) + "].MaxRecursionLevel", 3);
+ gmsh::option::setNumber("View[" + std::to_string(viewIndex) + "].TargetError", -0.0001);
+}
\ No newline at end of file
diff --git a/srcs/BEM/solverBEM.cpp b/srcs/BEM/solverBEM.cpp
index 248d0871dcc90e5c82737606429d3826e9bea724..67ab21db5341f66df856edbd53451b30359f963e 100644
--- a/srcs/BEM/solverBEM.cpp
+++ b/srcs/BEM/solverBEM.cpp
@@ -11,44 +11,29 @@
#include <omp.h>
#include <list>
-void solverBEM(std::map<int, double> &electrostaticPressure, int & nbViews, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, bool postProcessing, const int iteration, bool untangleMesh)
+// Fills the "electrostaticPressure" map that gives the electrostatic pressure associated with each (boundary) element tag.
+// "nbViews" is the counter of views (useful if other views have to be generated than with this solver). "nodalDisplacements" is
+// a map associating the (x,y) displacement to a serie of node tags (if the position of the nodes are not striclty the one
+// extracted from the .geo file). "postProcessing" indicates whether the post-processing should be done ('true') or not
+// ('false'). "iteration" represents the number of times the solver has been called (some information is printed on the first
+// iteration). "untangleMesh" unstangles every BEM domain if set to 'true' (in case the "nodalDisplacements" map contains
+// non-zero displacements).
+void solverBEM(std::map<int, double> &electrostaticPressure,
+ int & nbViews,
+ std::map<int,std::pair<double,double>> & nodalDisplacements,
+ const bool postProcessing,
+ const int iteration,
+ const bool untangleMesh)
{
+ bool computePhi = true; // true: computes the potential.
+ bool computeElectricField = true; // true: computes the electric field.
- if(untangleMesh)
- {
- /*--------get physical groups--------------*/
- std::map<std::string, std::pair<int, int>> groups;
- gmsh::vectorpair dimTags;
- gmsh::model::getPhysicalGroups(dimTags);
- for (std::size_t i = 0; i < dimTags.size(); ++i)
- {
- int dim = dimTags[i].first;
- int tag = dimTags[i].second;
- std::string name;
- gmsh::model::getPhysicalName(dim, tag, name);
- groups[name] = {dim, tag};
- }
-
- std::string groupname = "BEM_domain_1";
- int dim = groups[groupname].first;
- int tag = groups[groupname].second;
-
- std::vector<int> tags;
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
-
- gmsh::vectorpair dimTagsBis(tags.size());
- for(size_t i = 0; i < tags.size(); i++)
- {
- dimTagsBis[i] = std::pair<int,int>(dim, tags[i]);
- }
-
- gmsh::model::mesh::optimize("UntangleMeshGeometry", false, 1, dimTagsBis);
- }
-
+ // Retrieves the names of the BEM domains as well as their physical dim/tag.
gmsh::vectorpair physicalGroupsDimTags;
gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
- int nbDomains = 0;
+ std::vector<std::string> domainNames;
+ std::vector<std::pair<int, int>> domainPhysicalGroupDimTags;
for(std::size_t i = 0; i < physicalGroupsDimTags.size(); i++)
{
@@ -63,385 +48,247 @@ void solverBEM(std::map<int, double> &electrostaticPressure, int & nbViews, std:
while(getline(ss, word, '_'))
words.push_back(word);
- if(words.size() == 3 && words[0].compare("BEM") == 0)
+ if(words.size() == 3 && words[0].compare("BEM") == 0 && words[1].compare("domain") == 0)
{
- if(words[1].compare("domain") == 0)
- {
- if(stoi(words[2]) > nbDomains)
- {
- nbDomains = stoi(words[2]);
- }
- }
+ domainNames.push_back(words[2]);
+ domainPhysicalGroupDimTags.push_back({physicalGroupsDim, physicalGroupsTag});
}
}
- if(iteration == 1)
+ // Useful for dealing with the untangler.
+ if(untangleMesh)
+ // Untangles each BEM domain.
+ for(std::size_t i = 0; i < domainNames.size(); i++)
+ {
+ int dim = domainPhysicalGroupDimTags[i].first;
+ int tag = domainPhysicalGroupDimTags[i].second;
+
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+ gmsh::vectorpair dimTagsBis(tags.size());
+ for(size_t i = 0; i < tags.size(); i++)
+ dimTagsBis[i] = std::pair<int,int>(dim, tags[i]);
+
+ gmsh::model::mesh::optimize("UntangleMeshGeometry", false, 1, dimTagsBis);
+
+ }
+
+
+ if (domainNames.size() == 0)
{
- std::cout << "The model contains " << nbDomains << " BEM domain(s) \n";
+ std::cerr << "Error: No BEM domain found ! Usage: <BEM_domain_'domainName'>\n";
+ exit(0);
}
- if (nbDomains == 0)
+ if(iteration == 1)
{
- std::cerr << "Wrong names for the BEM domains. Usage: <BEM_domain_'i'>\n";
- exit(0);
+ std::cout << "--------------------[BEM SOLVER]--------------------\n";
+ std::cout << "The model contains " << domainNames.size() << " BEM domain(s):\n";
}
- //#pragma omp parallel for
- for(int i = 1; i <= nbDomains; i++)
+ // Generates viewTags only if the postProcessing is enabled.
+ int phiViewTag = 0;
+ int electricFieldViewTag = 0;
+ if(postProcessing)
{
- singleDomainSolverBEM(electrostaticPressure, boundaryDisplacementMap, nbViews, postProcessing, iteration, i);
+ if(computePhi)
+ phiViewTag = gmsh::view::add("Phi");
+ if(computeElectricField)
+ electricFieldViewTag = gmsh::view::add("Electric Field");
}
-}
-void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int & nbViews, bool postProcessing, const int iteration, const int domainTag)
-{
- // Parameters for displaying variables (for debugging).
- bool dispNodeCoord = 0; // Coordinates of the nodes.
- bool dispHierarchy = 0; // Hierarchy (physical groups -> entities -> elements -> nodes).
- bool dispMatrices = 0; // Matrices A, B, c, b, x. /!\ May lead to visualisation problems for large matrices.
- bool dispValues = 0; // Elements and their values of u and u_n.
- bool show_element_node_data = 0; //show the nodes and elements on the bondary that are in contact with the elements of the domain
- bool show_time = 0; //display the elapsed time for the internal potential calculus in the terminal
- bool show_node_class = 0; //shows the class of the node. Either if it belongs to the domain or to the boundary.
-
- // Code parameters to activate different kind of algorithms
- bool convergence = 0; // compute the convergence between the analytical solution and the numeric one /!\ use this routine with rectangle.geo
- bool electric_field = 1; //0: no electric field visualization, 1: compute and display the electric field in the domain
- int visualisation = 1; //0: no visualisation, 1: element node data, 2: element data
- bool solution_type = 1; //0: numeric, 1: analytique
-
- //to display the time of the different functions
- std::chrono::time_point<std::chrono::system_clock> start, end;
- std::chrono::time_point<std::chrono::system_clock> start_total, end_total;
- std::string function_name;
- start_total = std::chrono::system_clock::now();
-
- // Initialises a vector of elements (that will contain all the elements of interrest)
- std::vector<elementStruct> elementVector;
- std::string integrationType = "Gauss4"; // Could vary depending on the distance between the elements.
+ // Runs the solver for each BEM domain.
+ for(std::size_t i = 0; i < domainNames.size(); i++)
+ singleDomainSolverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, iteration, "BEM_domain_" + domainNames[i], phiViewTag, electricFieldViewTag, computePhi, computeElectricField);
- // Gets the tags, the coordinates and the parametric coordinates of the nodes.
- std::map<size_t, int> allNodeIndex; // nodeTag -> position in allNodeTags and allCoord
- std::vector<size_t> allNodeTags;
- std::vector<double> allCoord;
- std::vector<double> allParametricCoord;
- gmsh::model::mesh::getNodes(allNodeTags, allCoord, allParametricCoord);
+ if(iteration == 1)
+ std::cout << "--------------------[BEM SOLVER]--------------------\n\n";
- for(size_t i = 0; i < allNodeTags.size(); i++)
+ // Updates the number of views.
+ if(postProcessing)
{
- allNodeIndex[allNodeTags[i]] = i;
+ if(computePhi)
+ nbViews++;
+ if(computeElectricField)
+ nbViews++;
}
+}
- // Print the coordinates of all nodes if dispNodeCoord != 0
- if(dispNodeCoord)
- {
- dispNodeCoordFun(allNodeTags, allCoord);
- }
- // Initially, these are 0 and are incremented each time a dirichlet/neumann-element is encountered.
- int dirichletCount = 0;
-
- //pour compute_internal_phi
- std::vector<size_t> NodeTags2d_domain;
- std::vector<double> NodeCoord2d_domain;
- // vector<double> NodeParametricCoord2d_domain;
- // model::mesh::getNodes(NodeTags2d_domain, NodeCoord2d_domain, NodeParametricCoord2d_domain,2);
-
- //pour elementnodedata
- std::vector<int> elementTypes2d_domain_BEM;
- std::vector<std::vector<size_t> > elementTags2d_domain_BEM;
- std::vector<std::vector<size_t> > nodeTags2d_domain_BEM;
- int entities_tags_for_physical_group_domain_BEM;
- std::vector<int> entities_tags_2d_domain_BEM;
- std::map<std::size_t, std::pair<bool,bool>> nodeOnBoundary;
-
- // model::mesh::getElements(elementTypes2d,elementTags2d,nodeTags2d,2);
-
- start = std::chrono::system_clock::now();
- //calculates the potential at the barycentre of each element of the domain
- getBC(allNodeIndex, allCoord, elementVector, dirichletCount, dispHierarchy, NodeTags2d_domain, NodeCoord2d_domain,elementTypes2d_domain_BEM,elementTags2d_domain_BEM,nodeTags2d_domain_BEM,entities_tags_for_physical_group_domain_BEM, entities_tags_2d_domain_BEM, boundaryDisplacementMap, nodeOnBoundary,domainTag);
- end = std::chrono::system_clock::now();
- function_name = "getBC";
- display_time(start, end, function_name, show_time, domainTag);
-
- if(show_node_class)
- {
- for(std::size_t i = 0; i < NodeTags2d_domain.size(); i++)
- {
- std::cout << "node Tag: " << NodeTags2d_domain[i] << " is on the domain(0) or on the boundary(1): " <<
- nodeOnBoundary[NodeTags2d_domain[i]].first << std::endl;
- }
- }
+// Fills the "electrostaticPressure" map. The "nodalDisplacements", "postProcessing" and "iteration" arguments are described in
+// the description of the function "solverBEM". "domainName" is used for extracting everything in the appropriate independant BEM
+// domain, "phiViewTag" and "electricFieldViewTag" are the view tag of the potential and of the electric field respectively.
+// "computePhi" and "computeElectricField" indicates whether the potential and/or the electric field should be computed ('true')
+// or not ('false').
+void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
+ const int nbViews,
+ std::map<int,std::pair<double,double>> & nodalDisplacements,
+ const bool postProcessing,
+ const int iteration,
+ const std::string domainName,
+ const int phiViewTag,
+ const int electricFieldViewTag,
+ const bool computePhi,
+ const bool computeElectricField)
+{
+ // Parameters:
+ bool showTime = false; // true: displays the different time measurements.
+
+ bool convergence = false; // true: computes the error between the analytical solution and the numerical one on a simple configuration (/!\ use this with the BEM_convergence.geo model).
+ bool visualisation = true; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of the potential (Remark: electric field is always displayed with 'ElementData' type).
+ bool solutionType = true; // true: analytical computations, false: numerical computations.
+
+ // Used for displaying the time of different sections.
+ // double start, end;
+ double startTotal, endTotal;
+ double start, end;
+ startTotal = omp_get_wtime();
+
+ // Initializes a vector of elements (that will contain all the elements of interrest)
+ std::vector<elementStruct> elementVector;
+ std::string integrationType = "Gauss10"; // Used for the numerical integrations when "solutionType" = 'false'.
- std::vector<std::string> names;
- gmsh::model::list(names);
- //std::cout << "the file contains " << names.size() << " model names\n";
+ // Fills the vector of boundary elements.
+ start = omp_get_wtime();
- // Declaration of matrices A, B and vector c (the matrix type is perhaps not the most efficient).
+ std::map<std::string, std::pair<int, int>> physicalGroups;
+ std::size_t nbDirichletElements = fillElementVector(elementVector, physicalGroups, nodalDisplacements, domainName);
+
+ end = omp_get_wtime();
+ displayTime(start, end, "'fillElementVector'", showTime);
+
+ if(iteration == 1)
+ std::cout << "\n- BEM domain '" << domainName << "': contains " << elementVector.size() << " boundary elements.\n";
+
+ // Constructs the linear system
+ start = omp_get_wtime();
+
+ // The matrix type is perhaps not the most efficient.
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(elementVector.size(), elementVector.size());
Eigen::MatrixXd B = Eigen::MatrixXd::Zero(elementVector.size(), elementVector.size());
Eigen::MatrixXd c = Eigen::MatrixXd::Zero(elementVector.size(), 1);
+ fillSystem(A, B, c, elementVector, nbDirichletElements, integrationType, solutionType);
- // Double loop over each element in order to compute the element h and g that will be given to the matrices A and B.
- start = std::chrono::system_clock::now();
- double h;
- double g;
- #pragma omp parallel for schedule(guided) shared(elementVector,A,B,c,integrationType,dirichletCount) private(h,g)
- for(int i = 0; i < (int) elementVector.size(); i++)
- {
- for(int j = 0; j < (int) elementVector.size(); j++)
- {
- // Will be computed and stored in A and B.
- h = 0;
- g = 0;
+ end = omp_get_wtime();
+ displayTime(start, end, "'fillSystem'", showTime);
- if(i == j) // Formula for diagonal elements.
- {
- h = -0.5;
- g = (elementVector[i].length/(2*M_PI)) * (log(elementVector[i].length/2)-1);
- }
- else // General case
- {
- h = computeH(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
- if(solution_type == 0)
- {
- g = computeG(elementVector[i].midX, elementVector[i].midY, elementVector[j], integrationType);
- }
- else
- {
- g = computeG_analytique(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
- }
- }
+ // Solves the linear system.
+ start = omp_get_wtime();
- // Storage of h and g at the right place (in matrices A and B).
- if(j < dirichletCount)
- {
- A(i,j) = g;
- B(i,j) = h;
- c(j) = elementVector[j].bcValue[0];
- }
- else
- {
- A(i,j) = -h;
- B(i,j) = -g;
- c(j) = elementVector[j].bcValue[1];
- }
- }
- }
- end = std::chrono::system_clock::now();
- function_name = "compute H and G";
- display_time(start, end, function_name, show_time, domainTag);
-
-
- //If we want to display the time it takes for eigen to solve the linear system
- start = std::chrono::system_clock::now();
Eigen::MatrixXd b = B*c;
- Eigen::MatrixXd x = A.fullPivLu().solve(b);
- end = std::chrono::system_clock::now();
- function_name = "solve linear system";
- display_time(start, end, function_name, show_time, domainTag);
+ Eigen::MatrixXd x = A.partialPivLu().solve(b);
- // Matrices verification
- if(dispMatrices)
- {
- dispMatricesFun(A, B, c, b, x);
- }
+ end = omp_get_wtime();
+ displayTime(start, end, "Solve linear system", showTime);
- // Assign each value of x to the right element.
- for(int i = 0; i < (int) elementVector.size(); i++)
- {
- if( i < dirichletCount)
- elementVector[i].bcValue[1] = x(i);
- else
- elementVector[i].bcValue[0] = x(i);
- }
+ // Computes the electrostatic pressure on each boundary element.
+ start = omp_get_wtime();
+ computeElectrostaticPressure(electrostaticPressure, elementVector, x, nbDirichletElements, domainName);
+ end = omp_get_wtime();
+ displayTime(start, end, "'computeElectrostaticPressure'", showTime);
- if(dispValues)
+ // Visualization
+ if(postProcessing)
{
- dispValuesFun(elementVector);
- }
+ int viewIndex = nbViews;
+
+ // Retrieves everything related to domain of interrest.
+ start = omp_get_wtime();
- /*
- * Visualisation
- */
+ std::vector<std::size_t> domainNodeTags;
+ std::map<std::size_t, std::size_t> domainNodeIndices;
+ std::vector<double> domainNodeCoords;
- //visualisation by elementnodedata
- if(visualisation == 1 && postProcessing)
- {
- //The first index indicates the elementType, the second the elementTag and the third the nodesTags associated with this element. For example, it says whether it is a square or a triangle, the element tag and finally the tag of the 3 or 4 nodes associated with it.
- std::vector<std::vector<std::vector<double>>> data;
- //This map associates the value of the potential to all nodes strictly included in the domain and not on the boundary.
- std::map<int, double> nodeMapDomain;
-
- //To displays the time of the this code section
- //the purpose of this function is to fill the nodemapdomain map and thus to calculate the value of the potential at each node in the domain
- start = std::chrono::system_clock::now();
- //calculates the potential at the barycentre of each element of the domain
- compute_internal_phi(allCoord,NodeTags2d_domain,allNodeIndex,elementVector,nodeMapDomain,integrationType,nodeOnBoundary,solution_type);
- end = std::chrono::system_clock::now();
- function_name = "compute_internal_phi";
- display_time(start, end, function_name, show_time, domainTag);
-
- //this function allows to correctly associate the value of the potential to each node of an element in order to display the ElementNodeData view
- start = std::chrono::system_clock::now();
- element_node_data(data,elementTags2d_domain_BEM,elementTypes2d_domain_BEM,nodeTags2d_domain_BEM,elementVector, nodeMapDomain,show_element_node_data,entities_tags_for_physical_group_domain_BEM,nodeOnBoundary);
- end = std::chrono::system_clock::now();
- std::string function_name = "ElementNodeData";
- display_time(start, end, function_name, show_time, domainTag);
-
- #pragma omp critical (addPhiView)
- {
- int viewNodeTag = gmsh::view::add("phi");
+ std::vector<int> domainEntityTags;
+ std::map<int, std::size_t> domainEntityIndices;
+ std::vector<std::vector<int>> domainElementTypes;
+ std::vector<std::vector<std::vector<size_t>>> domainElementTags;
+ std::map<std::size_t, std::vector<std::size_t>> domainElementNodeTags;
- for(size_t i = 0; i < data.size(); i++)
- {
- gmsh::view::addModelData(viewNodeTag,0,names[0],"ElementNodeData",elementTags2d_domain_BEM[i],data[i],0.0,1);
- }
- gmsh::view::write(viewNodeTag, "results.msh");
-
- //gmsh::option::setNumber("View[" + std::to_string(domainTag-1) + "].IntervalsType", 3);
- //option::setNumber("View[0].AdaptVisualizationGrid", 1);
- //gmsh::option::setNumber("View[0].MaxRecursionLevel", 3);
- //gmsh::option::setNumber("View[0].TargetError", -0.0001);
- std::string my_string = "View[" + std::to_string(nbViews) + "].Visible";
- gmsh::option::setNumber(my_string, 0);
- nbViews++; // number of views generated (useful for FEM display of results)
- }
+ int domainPhysicalGroupTag = getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainEntityIndices, domainElementTypes, domainElementTags, domainElementNodeTags, physicalGroups, domainName);
+
+ end = omp_get_wtime();
+ displayTime(start, end, "'getModel'", showTime);
+
+ // The first index indicates the "entityTag", the second the "elementType", the third the "elementTag", and the fourth,
+ // either the nodesTags associated with the element (for elementNodalPhi) or the computed values (for "elemental...") (1
+ // for the potential and 3 for the electric field).
+ std::vector<std::vector<std::vector<std::vector<double>>>> elementalPhi;
+ std::vector<std::vector<std::vector<std::vector<double>>>> elementNodalPhi;
+ std::vector<std::vector<std::vector<std::vector<double>>>> elementalElectricField;
- if(electric_field)
- {
- int nb2dElements = 0;
- for(size_t i = 0; i < elementTypes2d_domain_BEM.size(); i++)
- {
- nb2dElements += elementTags2d_domain_BEM[i].size();
- }
- std::vector<size_t> tags(nb2dElements);
-
- start = std::chrono::system_clock::now();
- std::vector<std::vector<double>> data_Electric_field(nb2dElements,std::vector<double>(3,0));
- electricField(nb2dElements,entities_tags_2d_domain_BEM,elementVector,tags,integrationType, data_Electric_field,solution_type);
- end = std::chrono::system_clock::now();
- function_name = "compute Electric Field";
- display_time(start, end, function_name, show_time, domainTag);
- #pragma omp critical (addElecView)
- {
- int viewtag = gmsh::view::add("Electric field");
- gmsh::view::addModelData(viewtag, 0, names[0], "ElementData",tags, data_Electric_field);
- //gmsh::view::write(viewtag, "results.msh");
- //gmsh::option::setNumber("View[0].IntervalsType", 3);
- //gmsh::option::setNumber("View[0].MaxRecursionLevel", 3);
- //gmsh::option::setNumber("View[0].TargetError", -0.0001);
- std::string my_string = "View[" + std::to_string(nbViews) + "].Visible";
- gmsh::option::setNumber(my_string, 0);
- nbViews++;
- }
- }
- }
- //visualisation by element barycenter
- if(visualisation == 2 && postProcessing)
- {
std::vector<std::string> names;
gmsh::model::list(names);
- //To get the total number of 2d elements that make up the domain, we simply add the number of squares to the number of triangles,...
- int nb2dElements = 0;
- for(size_t i = 0; i < elementTypes2d_domain_BEM.size(); i++)
- nb2dElements += elementTags2d_domain_BEM[i].size();
-
- //Tag and associated potential of all 2d elements. In this view the potential is calculated at the barycentre and is constant per element.
- std::vector<std::vector<double>> phi(nb2dElements,std::vector<double>(1,0));
- std::vector<size_t> tags(nb2dElements);
-
- //To displays the time of the this code section
- start = std::chrono::system_clock::now();
- //calculates the potential at the barycentre of each element of the domain
- computeData(elementVector, tags, phi, integrationType, entities_tags_2d_domain_BEM,solution_type);
- end = std::chrono::system_clock::now();
- function_name = "computeData";
- display_time(start, end, function_name, show_time, domainTag);
-
- if(electric_field)
+ // Visualisation by ElementNodeData
+ if(visualisation)
{
- start = std::chrono::system_clock::now();
- std::vector<std::vector<double>> data_Electric_field(nb2dElements,std::vector<double>(3,0));
- electricField(nb2dElements,entities_tags_2d_domain_BEM,elementVector,tags,integrationType, data_Electric_field,solution_type);
- end = std::chrono::system_clock::now();
- function_name = "compute Electric Field";
- display_time(start, end, function_name, show_time, domainTag);
- #pragma omp critical (addElecView)
+ // Computes a map that stores 'nodeTag' (of the boundary) -> vector of adjacent 'elementIndices' (because only 1 element for higher order elements).
+ std::map<std::size_t, std::vector<std::size_t>> boundaryNodeToElements;
+ computeNodeToElements(boundaryNodeToElements, elementVector);
+
+ // Computes the value of the potential at each node strictly included in the domain and stores it in the "internalNodalPhi" map
+ start = omp_get_wtime();
+
+ std::map<int, double> internalNodalPhi;
+ computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector, integrationType, boundaryNodeToElements, solutionType);
+
+ end = omp_get_wtime();
+ displayTime(start, end, "'computeInternalPhi'", showTime);
+
+ // Associates the value of the potential to each node of an element.
+ start = omp_get_wtime();
+ elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector, internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements);
+ end = omp_get_wtime();
+ displayTime(start, end, "'elementNodeData'", showTime);
+
+ // Computes the electric field (with "ElementData" type).
+ start = omp_get_wtime();
+ computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, false, computeElectricField);
+ end = omp_get_wtime();
+ displayTime(start, end, "'computeElementData' (only the electric field)", showTime);
+
+ // Displays what is needed.
+ if(computePhi)
{
- int viewtag = gmsh::view::add("Electric field");
- gmsh::view::addModelData(viewtag, 0, names[0], "ElementData",tags, data_Electric_field);
- //gmsh::view::write(viewtag, "results.msh");
- //gmsh::option::setNumber("View[0].IntervalsType", 3);
- //gmsh::option::setNumber("View[0].MaxRecursionLevel", 3);
- //gmsh::option::setNumber("View[0].TargetError", -0.0001);
- // empeche que tout se plot l'un sur l'autre
- std::string my_string = "View[" + std::to_string(nbViews) + "].Visible";
- gmsh::option::setNumber(my_string, 0);
-
- nbViews++;
+ displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi, viewIndex);
+ viewIndex++;
}
+ if(computeElectricField)
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex);
}
-
- /*
- * Convergence
- */
-
- if(convergence)
- {
- convergenceFun(entities_tags_2d_domain_BEM,tags,phi);
- }
-
- //For ease of use, a function handles the code for the display.
- #pragma omp critical (addPhiView)
- {
- displayData("phi", names[0], "results.msh", tags, phi);
- std::string my_string = "View[" + std::to_string(nbViews) + "].Visible";
- gmsh::option::setNumber(my_string, 0);
- nbViews++; // number of views generated (useful for FEM display of results)
- }
- }
-
- // retrieving the dielectric permittivity of the particular BEM domain
- std::vector<std::string> keys;
- gmsh::onelab::getNames(keys, "(Materials).+");
- std::string domainName = "BEM_domain_" + std::to_string(domainTag);
- double epsilon = 0; // dielectric permittivity
- for (auto &key : keys)
- {
- // get corresponding value
- std::vector<double> value;
- gmsh::onelab::getNumber(key, value);
- // expected key structure is "type/group_name/field"
- // => split the key string into components
- std::stringstream ss(key);
- std::vector<std::string> words;
- std::string word;
- while(std::getline(ss, word, '/')) // read string until '/'
- words.push_back(word);
- if(words.size()==3)
+ // Visualisation by ElementData
+ else
{
- if(words[2].compare("Epsilon") == 0 && words[1].compare(domainName) == 0){
- epsilon = value[0];
+ // Computes the potential and the electric field.
+ start = omp_get_wtime();
+ computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, computePhi, computeElectricField);
+ end = omp_get_wtime();
+ displayTime(start, end, "'computeElementData'", showTime);
+
+ // Displays what is needed.
+ if(computePhi)
+ {
+ displayData(phiViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalPhi, viewIndex);
+ viewIndex++;
+ }
+ if(computeElectricField)
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex);
+
+ // Study of the error on a precise model.
+ if(convergence)
+ {
+ std::cout << "\tCONVERGENCE:\n";
+ std::string convergenceModelName = "BEM_convergence";
+ if(names[0].compare(convergenceModelName) == 0)
+ convergenceFun(domainEntityTags, domainElementTypes, domainElementTags, elementalPhi, elementalElectricField);
+ else
+ std::cout << "\tWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '" << convergenceModelName + ".geo" << "' model if you want to check for convergence.\n\n";
}
}
}
- //std::cout << "epsilon test: " << epsilon << "\n";
- if(epsilon == 0) // quid si solveur n'est pas couplé ?
- {
- std::cout << "you forgot to set the dielectric permittivity in the .geo file. \n";
- }
-
- //Calculates the electrostatic pressure per element in order to provide these results to the FEM code that will handle the rest of the algorithm
- //double epsilon = 8.8541878128e-12; // Ã tirer du .geo
- for(size_t i = 0; i < elementVector.size(); i++)
- {
- #pragma omp atomic write
- electrostaticPressure[elementVector[i].tag] = epsilon * pow(elementVector[i].bcValue[1], 2.0) / 2; // not forget factor 2
- }
- end_total = std::chrono::system_clock::now();
- function_name = "total CPU time";
- display_time(start_total, end_total, function_name, show_time, domainTag);
+ endTotal = omp_get_wtime();
+ displayTime(startTotal, endTotal, "total time of the code BEM", showTime);
}
\ No newline at end of file
diff --git a/srcs/BEM/verification.cpp b/srcs/BEM/verification.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48d7bbde1138e9dd2ccb7bd97b1502ea3555a34e
--- /dev/null
+++ b/srcs/BEM/verification.cpp
@@ -0,0 +1,109 @@
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+#include <sstream>
+#include <Eigen/Dense>
+#include <cmath>
+#include <iomanip>
+#include <stdio.h>
+#include <chrono>
+#include "functionsBEM.hpp"
+#include <omp.h>
+#include <algorithm>
+#include <list>
+
+void displayTime(const double start,
+ const double end,
+ std::string functionName, const bool showTime)
+{
+ //displays the time of the desired code section
+ if(showTime)
+ {
+ double duration = end - start;
+ std::cout << functionName << "\n";
+ std::cout << "Elapsed time: " << 1000000*duration << " [µs]\tor " << duration << " [s]\n\n";
+
+ // std::cout << std::endl;
+ // std::cout << function_name << std::endl;
+ // std::cout << "Elapsed time in nanoseconds: "
+ // << std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count()
+ // << " ns" << std::endl;
+
+ // std::cout << "Elapsed time in microseconds: "
+ // << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
+ // << " µs" << std::endl;
+
+ // std::cout << "Elapsed time in milliseconds: "
+ // << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
+ // << " ms" << std::endl;
+
+ // std::cout << "Elapsed time in seconds: "
+ // << std::chrono::duration_cast<std::chrono::seconds>(end - start).count()
+ // << " sec";
+ // std::cout << std::endl;
+ // std::cout << std::endl;
+ }
+}
+
+void convergenceFun(const std::vector<int> &domainEntityTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalPhi,
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalElectricField)
+{
+ double l = 1.0;
+ double phiLeft = 100.0;
+ double phiRight = 200.0;
+ double exactE = (phiLeft - phiRight) / l;
+
+ double xLeft = 0.25;
+ double xRight = 0.75;
+ double yBottom = 0.25;
+ double yTop = 0.75;
+
+ double errorPhi = 0.0;
+ double errorE = 0.0;
+ double errorEWindow = 0.0;
+ int count = 0;
+ int countWindow = 0;
+
+ std::vector<double> barycenters;
+
+ for(std::size_t i = 0; i < domainEntityTags.size(); i++)
+ {
+ for(std::size_t j = 0; j < domainElementTypes[i].size(); j++)
+ {
+ gmsh::model::mesh::getBarycenters(domainElementTypes[i][j], domainEntityTags[i], false, false, barycenters);
+ for(std::size_t k = 0; k < domainElementTags[i][j].size(); k++)
+ {
+ double numericalPhi = elementalPhi[i][j][k][0];
+ double numericalE = elementalElectricField[i][j][k][0];
+
+ double x = barycenters[3*k];
+ double y = barycenters[3*k + 1];
+
+ double exactPhi = phiLeft + (phiRight - phiLeft) * x / l;
+
+ double relErrorPhi = abs( (exactPhi - numericalPhi) / exactPhi );
+ double relErrorE = abs( (exactE - numericalE) / exactE);
+
+ if(x > xLeft && x < xRight && y > yBottom && y < yTop)
+ {
+ errorEWindow += relErrorE*relErrorE;
+ countWindow++;
+ }
+
+ errorPhi += relErrorPhi*relErrorPhi;
+ errorE += relErrorE*relErrorE;
+ count++;
+ }
+ }
+ }
+ errorPhi = sqrt(errorPhi / count);
+ errorE = sqrt(errorE / count);
+ errorEWindow = sqrt(errorEWindow / countWindow);
+
+ std::cout << "\tRelative error on the electric potential: " << errorPhi << ".\n";
+ std::cout << "\tRelative error on the electric field: " << errorE << ".\n";
+ std::cout << "\tRelative error on the electric field (without the edges of the domain): " << errorEWindow << ".\n\n";
+}
\ No newline at end of file
diff --git a/srcs/CMakeLists.txt b/srcs/CMakeLists.txt
index 316d105f99d1ebcf39d949786063fd681853045a..89f6178972ca147bcfbb9ea53a40562d82f2b105 100644
--- a/srcs/CMakeLists.txt
+++ b/srcs/CMakeLists.txt
@@ -54,7 +54,7 @@ ENDIF()
# ------------------------------------------------------------------------------
-FILE(GLOB SRCS *.cpp FEM/solverFEM.cpp FEM/functionsFEM.cpp BEM/solverBEM.cpp BEM/functionsBEM.cpp)
+FILE(GLOB SRCS *.cpp FEM/solverFEM.cpp FEM/functionsFEM.cpp BEM/solverBEM.cpp BEM/boundaryElements.cpp BEM/postProcessing.cpp BEM/verification.cpp BEM/computations.cpp)
ADD_EXECUTABLE(solver ${SRCS})
TARGET_LINK_LIBRARIES(solver PUBLIC ${GMSH_LIBRARIES})
diff --git a/srcs/beamActuation.geo b/srcs/beamActuation.geo
new file mode 100644
index 0000000000000000000000000000000000000000..a37ee6d35de807819cba90a988a91c0491262946
--- /dev/null
+++ b/srcs/beamActuation.geo
@@ -0,0 +1,177 @@
+scale = 1e-6;
+
+bt = 1; // Beam tickness
+bh = 1; // Beam height (position)
+
+eL = 3; // Electrode base length
+el = 2.5; // Electrode tip length
+et = 1; // Electrode tickness
+es = 1.5; // Electrode spacing
+
+W = 3*es + 3*eL; // Width of the space
+H = bh + 3*(bt+et); // Height of the space
+
+nEs = 9;
+nEL = 12;
+nBt = 6;
+nEt = 3;
+nBh = 3;
+
+bt = bt * scale;
+bh = bh * scale;
+eL = eL * scale;
+el = el * scale;
+et = et * scale;
+es = es * scale;
+W = W * scale;
+H = H * scale;
+
+e = (eL - el)/2; // Electrode constant
+
+/*
+ * Points
+ */
+
+// Space
+Point(1) = {0, 0, 0, 1};
+Point(2) = {W, 0, 0, 1};
+Point(3) = {W, H, 0, 1};
+Point(4) = {0, H, 0, 1};
+
+// Beam
+Point(5) = {0, bh+bt, 0, 1};
+Point(6) = {es, bh+bt, 0, 1};
+Point(7) = {es+eL, bh+bt, 0, 1};
+Point(8) = {2*es+eL, bh+bt, 0, 1};
+Point(9) = {2*es+2*eL, bh+bt, 0, 1};
+Point(10) = {3*es+2*eL, bh+bt, 0, 1};
+
+Point(11) = {3*es+2*eL, bh, 0, 1};
+Point(12) = {2*es+2*eL, bh, 0, 1};
+Point(13) = {2*es+eL, bh, 0, 1};
+Point(14) = {es+eL, bh, 0, 1};
+Point(15) = {es, bh, 0, 1};
+Point(16) = {0, bh, 0, 1};
+
+// Electrodes
+Point(17) = {es+e, bh+bt+et, 0, 1};
+Point(18) = {es+e+el, bh+bt+et, 0, 1};
+Point(19) = {2*es+eL+e, bh+bt+et, 0, 1};
+Point(20) = {2*es+eL+e+el, bh+bt+et, 0, 1};
+
+/*
+ * Lines
+ */
+
+ // Space
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 5};
+Line(5) = {16, 1};
+
+// Beam
+Line(6) = {5, 6};
+Line(7) = {6, 7};
+Line(8) = {7, 8};
+Line(9) = {8, 9};
+Line(10) = {9, 10};
+Line(11) = {10, 11};
+Line(12) = {11, 12};
+Line(13) = {12, 13};
+Line(14) = {13, 14};
+Line(15) = {14, 15};
+Line(16) = {15, 16};
+Line(17) = {16, 5};
+Line(18) = {15, 6};
+Line(19) = {14, 7};
+Line(20) = {13, 8};
+Line(21) = {12, 9};
+
+// Electrode
+Line(22) = {6, 17};
+Line(23) = {17, 18};
+Line(24) = {18, 7};
+Line(25) = {8, 19};
+Line(26) = {19, 20};
+Line(27) = {20, 9};
+
+/*
+ * Transfinite Curve
+ */
+Transfinite Curve{6, 8, 10, 12, 14, 16} = nEs + 1 Using Progression 1;
+Transfinite Curve{7, 9, 13, 15, 23, 26} = nEL + 1 Using Progression 1;
+Transfinite Curve{17, 18, 19, 20, 21, 11} = nBt + 1 Using Progression 1;
+Transfinite Curve{22, 24, 25, 27} = nEt + 1 Using Progression 1;
+
+Transfinite Curve{1, 3} = 3*nEs + 3*nEL + 1 Using Progression 1;
+Transfinite Curve(2) = 4 * (nBt + nEt) + 1 Using Progression 1;
+Transfinite Curve(4) = 3*nEt + 2*nBt + 1 Using Progression 1;
+Transfinite Curve(5) = nBh + 1 Using Progression 1;
+
+/*
+ * Curve Loop
+ */
+Curve Loop(1) = {6, -18, 16, 17};
+Curve Loop(2) = {7, -19, 15, 18};
+Curve Loop(3) = {8, -20, 14, 19};
+Curve Loop(4) = {9, -21, 13, 20};
+Curve Loop(5) = {10, 11, 12, 21};
+Curve Loop(6) = {23, 24, -7, 22};
+Curve Loop(7) = {26, 27, -9, 25};
+Curve Loop(8) = {1, 2, 3, 4, 6, 22, 23, 24, 8, 25, 26, 27, 10, 11, 12, 13, 14, 15, 16, 5};
+
+/*
+ * Plane Surface
+ */
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {3};
+Plane Surface(4) = {4};
+Plane Surface(5) = {5};
+Plane Surface(6) = {6};
+Plane Surface(7) = {7};
+Plane Surface(8) = {8};
+
+/*
+ * Surface Manipulation
+ */
+Transfinite Surface{1, 2, 3, 4, 5, 6, 7};
+Recombine Surface{1, 2, 3, 4, 5, 6, 7};
+
+/*
+ * Physical Curve
+ */
+Physical Curve("clamp", 1) = {17};
+Physical Curve("outside", 2) = {1, 2, 3, 4, 5};
+Physical Curve("left_electrode", 3) = {24};
+Physical Curve("right_electrode", 4) = {25};
+Physical Curve("beam", 5) = {6, 22, 23, 8, 26, 10, 11, 12, 13, 14, 15, 16};
+Physical Curve("BEM_FEM_boundary", 6) = {6, 22, 23, 24, 8, 26, 26, 27, 10, 11, 12, 13, 14, 15, 16};
+
+/*
+ * Physical Surface
+ */
+Physical Surface("FEM_domain", 1) = {1, 2, 3, 4, 5, 6, 7};
+Physical Surface("BEM_domain_1", 2) = {8};
+
+/*
+ * Parameters
+ */
+SetNumber("Boundary Conditions/clamp/ux", 0.0);
+SetNumber("Boundary Conditions/clamp/uy", 0.0);
+
+SetNumber("Materials/FEM_domain/Young", 3e7);
+SetNumber("Materials/FEM_domain/Poisson", 0.3);
+SetNumber("Materials/FEM_domain/rho",7800);
+
+SetNumber("Volumic Forces/FEM_domain/bx", 0.0);
+SetNumber("Volumic Forces/FEM_domain/by", 0.0);
+
+SetNumber("Boundary Conditions/outside/BEM_domain_1/neumann", 0);
+SetNumber("Boundary Conditions/beam/BEM_domain_1/neumann", 0);
+SetNumber("Boundary Conditions/left_electrode/BEM_domain_1/dirichlet", 200);
+SetNumber("Boundary Conditions/right_electrode/BEM_domain_1/dirichlet", 0);
+SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
+
+SetNumber("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/beamActuationGeneral.geo b/srcs/beamActuationGeneral.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3ad6cf8a80d55447284b1467115cedc983abd0e0
--- /dev/null
+++ b/srcs/beamActuationGeneral.geo
@@ -0,0 +1,229 @@
+scale = 1e-6;
+
+bt = 0.2; // Beam tickness
+bh = 1; // Beam height (position)
+
+eL = 3; // Electrode base length
+el = 1.5; // Electrode tip length
+et = 1; // Electrode tickness
+es = 1.5; // Electrode spacing
+
+nbElectrodes = 10;
+
+W = nbElectrodes*(es + eL) + 2*es; // Width of the space
+H = W; // Height of the space
+
+nEs = 12;
+nEL = 15;
+nBt = 8;
+nEt = 3;
+nBh = 3;
+
+bt = bt * scale;
+bh = bh * scale;
+eL = eL * scale;
+el = el * scale;
+et = et * scale;
+es = es * scale;
+W = W * scale;
+H = H * scale;
+
+e = (eL - el)/2; // Electrode constant
+
+/*
+ * Points
+ */
+
+For i In {0:nbElectrodes-1}
+ Point(6*i + 1) = {(es+eL)*i, bh + bt, 0, 1};
+ Point(6*i + 2) = {(es+eL)*i + es, bh + bt, 0, 1};
+ Point(6*i + 3) = {(es+eL)*i + es+e, bh + bt + et, 0, 1};
+ Point(6*i + 4) = {(es+eL)*i + es+e+el, bh + bt + et, 0, 1};
+ Point(6*i + 5) = {(es+eL)*i + es, bh, 0, 1};
+ Point(6*i + 6) = {(es+eL)*i, bh, 0, 1};
+EndFor
+Point(6*nbElectrodes + 1) = {nbElectrodes*(es+eL), bh + bt, 0, 1};
+Point(6*nbElectrodes + 2) = {nbElectrodes*(es+eL) + es, bh + bt, 0, 1};
+Point(6*nbElectrodes + 3) = {nbElectrodes*(es+eL) + es, bh, 0, 1};
+Point(6*nbElectrodes + 6) = {nbElectrodes*(es+eL), bh, 0, 1};
+
+Point(6*nbElectrodes + 4) = {0, 0, 0, 1};
+Point(6*nbElectrodes + 5) = {W, 0, 0, 1};
+Point(6*nbElectrodes + 7) = {W, H, 0, 1};
+Point(6*nbElectrodes + 8) = {0, H, 0, 1};
+
+For i In {0:nbElectrodes-1}
+ Point(6*nbElectrodes + 8 + i + 1) = {i*(es+eL) + es + e + el/2, bh + bt + et - ((el*e)/(2*et)), 0, 1};
+EndFor
+
+/*
+ * Lines
+ */
+Line(9*nbElectrodes + 1) = {6*nbElectrodes + 4, 6*nbElectrodes + 5};
+Line(9*nbElectrodes + 2) = {6*nbElectrodes + 5, 6*nbElectrodes + 7};
+Line(9*nbElectrodes + 3) = {6*nbElectrodes + 7, 6*nbElectrodes + 8};
+Line(9*nbElectrodes + 4) = {6*nbElectrodes + 8, 1};
+Line(9*nbElectrodes + 5) = {6, 6*nbElectrodes + 4};
+
+Line(9*nbElectrodes + 9) = {6*nbElectrodes + 1, 6*nbElectrodes + 2};
+Line(9*nbElectrodes + 6) = {6*nbElectrodes + 2, 6*nbElectrodes + 3};
+Line(9*nbElectrodes + 7) = {6*nbElectrodes + 3, 6*nbElectrodes + 6};
+Line(9*nbElectrodes + 8) = {6*nbElectrodes + 6, 6*nbElectrodes + 1};
+For i In {0:nbElectrodes-1}
+ Line(9*i + 1) = {6*i + 1, 6*i + 2};
+ Line(9*i + 2) = {6*i + 2, 6*i + 3};
+ // Line(9*i + 3) = {6*i + 3, 6*i + 4};
+ Circle(9*i + 3) = {6*i + 3, 6*nbElectrodes + 8 + i + 1, 6*i + 4};
+ Line(9*i + 4) = {6*i + 4, 6*(i+1) + 1};
+ Line(9*i + 5) = {6*(i+1) + 1, 6*i + 2};
+ Line(9*i + 6) = {6*i + 5, 6*i + 2};
+ Line(9*i + 7) = {6*i + 5, 6*i + 6};
+ Line(9*i + 8) = {6*i + 6, 6*i + 1};
+ Line(9*i + 9) = {6*(i+1) + 6, 6*i + 5};
+EndFor
+
+/*
+ * Transfinite Curve
+ */
+For i In {0:nbElectrodes-1}
+ Transfinite Curve{9*i + 1, 9*i + 7} = nEs + 1 Using Progression 1;
+ Transfinite Curve{9*i + 2, 9*i + 4} = nEt + 1 Using Progression 1;
+ Transfinite Curve{9*i + 3, 9*i + 5, 9*i + 9} = nEL + 1 Using Progression 1;
+ Transfinite Curve{9*i + 8, 9*i + 6} = nBt + 1 Using Progression 1;
+EndFor
+Transfinite Curve{9*nbElectrodes + 9, 9*nbElectrodes + 7} = nEs + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 8, 9*nbElectrodes + 6} = nBt + 1 Using Progression 1;
+
+Transfinite Curve{9*nbElectrodes + 1, 9*nbElectrodes + 3} = nbElectrodes*(nEs + nEL) + 2*nEs + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 2} = nbElectrodes*(nEs + nEL) + 2*nEs + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 4} = nbElectrodes*(nEs + nEL) + 2*nEs - nBt - nBh + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 5} = nBh + 1 Using Progression 1;
+
+/*
+ * Curve Loop
+ */
+
+For i In {0:nbElectrodes-1}
+ Curve Loop(3*i + 1) = {9*i + 1, -(9*i + 6), 9*i + 7, 9*i + 8};
+ Curve Loop(3*i + 2) = {9*i + 3, 9*i + 4, 9*i + 5, 9*i + 2};
+ Curve Loop(3*i + 3) = {-(9*i + 5), -(9*(i+1) + 8), 9*i + 9, 9*i + 6};
+EndFor
+Curve Loop(3*nbElectrodes + 1) = {9*nbElectrodes + 9, 9*nbElectrodes + 6, 9*nbElectrodes + 7, 9*nbElectrodes + 8};
+// How to fill a vector in a loop ?
+
+bemBoundary = {1:6*nbElectrodes + 8};
+bemBoundaryWithoutE = {1:4*nbElectrodes + 8 + 2};
+bemFemBoundary = {1:6*nbElectrodes + 3};
+leftElectrodesFull = {1:nbElectrodes};
+rightElectrodesFull = {1:nbElectrodes};
+For i In {0:nbElectrodes-1}
+ bemBoundary[4*i+1 - 1] = 9*i + 1;
+ bemBoundary[4*i+2 - 1] = 9*i + 2;
+ bemBoundary[4*i+3 - 1] = 9*i + 3;
+ bemBoundary[4*i+4 - 1] = 9*i + 4;
+
+ bemBoundaryWithoutE[2*i+1 - 1] = 9*i + 1;
+ bemBoundaryWithoutE[2*i+2 - 1] = 9*i + 3;
+
+ bemFemBoundary[4*i+1 - 1] = 9*i + 1;
+ bemFemBoundary[4*i+2 - 1] = 9*i + 2;
+ bemFemBoundary[4*i+3 - 1] = 9*i + 3;
+ bemFemBoundary[4*i+4 - 1] = 9*i + 4;
+
+ leftElectrodesFull[i] = 9*i + 2;
+ rightElectrodesFull[i] = 9*i + 4;
+EndFor
+bemBoundary[4*nbElectrodes + 1 - 1] = 9*nbElectrodes + 9;
+bemBoundary[4*nbElectrodes + 2 - 1] = 9*nbElectrodes + 6;
+bemBoundary[4*nbElectrodes + 3 - 1] = 9*nbElectrodes + 7;
+
+bemBoundaryWithoutE[2*nbElectrodes + 1 - 1] = 9*nbElectrodes + 9;
+bemBoundaryWithoutE[2*nbElectrodes + 2 - 1] = 9*nbElectrodes + 6;
+bemBoundaryWithoutE[2*nbElectrodes + 3 - 1] = 9*nbElectrodes + 7;
+
+bemFemBoundary[4*nbElectrodes + 1 - 1] = 9*nbElectrodes + 9;
+bemFemBoundary[4*nbElectrodes + 2 - 1] = 9*nbElectrodes + 6;
+bemFemBoundary[4*nbElectrodes + 3 - 1] = 9*nbElectrodes + 7;
+For i In {0:nbElectrodes-1}
+ bemBoundary[4*nbElectrodes + 3 + 2*i+1 - 1] = 9*i + 9;
+ bemBoundary[4*nbElectrodes + 3 + 2*i+2 - 1] = 9*i + 7;
+
+ bemBoundaryWithoutE[2*nbElectrodes + 3 + 2*i+1 - 1] = 9*i + 9;
+ bemBoundaryWithoutE[2*nbElectrodes + 3 + 2*i+2 - 1] = 9*i + 7;
+
+ bemFemBoundary[4*nbElectrodes + 3 + 2*i+1 - 1] = 9*i + 9;
+ bemFemBoundary[4*nbElectrodes + 3 + 2*i+2 - 1] = 9*i + 7;
+EndFor
+bemBoundary[6*nbElectrodes + 4 - 1] = 9*nbElectrodes + 5;
+bemBoundary[6*nbElectrodes + 5 - 1] = 9*nbElectrodes + 1;
+bemBoundary[6*nbElectrodes + 6 - 1] = 9*nbElectrodes + 2;
+bemBoundary[6*nbElectrodes + 7 - 1] = 9*nbElectrodes + 3;
+bemBoundary[6*nbElectrodes + 8 - 1] = 9*nbElectrodes + 4;
+
+bemBoundaryWithoutE[4*nbElectrodes + 4 - 1] = 9*nbElectrodes + 5;
+bemBoundaryWithoutE[4*nbElectrodes + 5 - 1] = 9*nbElectrodes + 1;
+bemBoundaryWithoutE[4*nbElectrodes + 6 - 1] = 9*nbElectrodes + 2;
+bemBoundaryWithoutE[4*nbElectrodes + 7 - 1] = 9*nbElectrodes + 3;
+bemBoundaryWithoutE[4*nbElectrodes + 8 - 1] = 9*nbElectrodes + 4;
+
+leftElectrodes = {1:nbElectrodes-1};
+rightElectrodes = {1:nbElectrodes-1};
+For i In {0:nbElectrodes-2}
+ leftElectrodes[i] = leftElectrodesFull[i+1];
+ rightElectrodes[i] = rightElectrodesFull[i];
+EndFor
+
+bemBoundaryWithoutE[4*nbElectrodes + 9 - 1] = 2;
+bemBoundaryWithoutE[4*nbElectrodes + 10 - 1] = 9*nbElectrodes + 4;
+
+Curve Loop(3*nbElectrodes + 2) = bemBoundary[];
+
+/*
+ * Plane Surface
+ */
+For i In {0:nbElectrodes-1}
+ Plane Surface(3*i + 1) = {3*i + 1};
+ Plane Surface(3*i + 2) = {3*i + 2};
+ Plane Surface(3*i + 3) = {3*i + 3};
+EndFor
+Plane Surface(3*nbElectrodes + 1) = {3*nbElectrodes + 1};
+Plane Surface(3*nbElectrodes + 2) = {3*nbElectrodes + 2};
+
+/*
+ * Surface operations
+ */
+Recombine Surface{1:3*nbElectrodes + 1};
+Transfinite Surface{1:3*nbElectrodes + 1};
+
+/*
+ * Physical Curve
+ */
+Physical Curve("clamp", 1) = {8};
+Physical Curve("BEM_boundary_without_e", 2) = bemBoundaryWithoutE[];
+Physical Curve("BEM_FEM_boundary", 3) = bemFemBoundary[];
+Physical Curve("left_electrodes", 4) = leftElectrodes[];
+Physical Curve("right_electrodes", 5) = rightElectrodes[];
+
+Physical Surface("FEM_domain") = {1:3*nbElectrodes + 1};
+Physical Surface("BEM_domain_1") = {3*nbElectrodes + 2};
+
+/*
+ * Parameters
+ */
+ SetNumber("Boundary Conditions/clamp/ux", 0.0);
+ SetNumber("Boundary Conditions/clamp/uy", 0.0);
+
+// SetNumber("Materials/FEM_domain/Young", 3e7);
+SetNumber("Materials/FEM_domain/Young", 1e7);
+ SetNumber("Materials/FEM_domain/Poisson", 0.3);
+ SetNumber("Materials/FEM_domain/rho",7800);
+
+ SetNumber("Volumic Forces/FEM_domain/bx", 0.0);
+ SetNumber("Volumic Forces/FEM_domain/by", 0.0);
+
+ SetNumber("Boundary Conditions/BEM_boundary_without_e/BEM_domain_1/neumann", 0);
+ SetNumber("Boundary Conditions/left_electrodes/BEM_domain_1/dirichlet", 8);
+ SetNumber("Boundary Conditions/right_electrodes/BEM_domain_1/dirichlet", 0);
+ SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
+
+ SetNumber("Non_linear_solver", 1);
\ No newline at end of file