diff --git a/srcs/BEM/boundaryElements.cpp b/srcs/BEM/boundaryElements.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7884bf59fbe7fa6c57c1b97b7fb2db77eff90e99
--- /dev/null
+++ b/srcs/BEM/boundaryElements.cpp
@@ -0,0 +1,346 @@
+#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 getBC(std::map<std::string, std::pair<bool, double>> &bcMap,
+ const std::string domainName)
+{
+ // Gets the boundary conditions from onelab and fills a map: "name of the physical group" -> (true/false, value).
+ // true: dirichlet boundary condition
+ // false: neumann boundary condition
+ 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 BEM boundary conditions.
+ if(words.size() == 4 && words[2].compare(domainName) == 0)
+ {
+ if(words[3].compare("dirichlet") == 0)
+ bcMap[words[1]] = {true, values[0]};
+ else if(words[3].compare("neumann") == 0)
+ bcMap[words[1]] = {false, values[0]};
+ }
+ }
+}
+
+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;
+
+ std::map<std::string, std::pair<bool, double>> bcMap;
+ getBC(bcMap, 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(bcMap.find(physicalGroupName) == 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" !
+ bool bcType = bcMap[physicalGroupName].first;
+ double bcValue = bcMap[physicalGroupName].second;
+
+ 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};
+ for(std::size_t j = 0; j < elementTypes.size(); j++)
+ {
+ int nbNodesPerElement = 2;
+ if(elementTypes[j] == 8)
+ nbNodesPerElement = 3;
+
+ std::vector<std::size_t> elementTags;
+ std::vector<std::size_t> elementNodeTags;
+ gmsh::model::mesh::getElementsByType(elementTypes[j], elementTags, elementNodeTags, entityTags[i]);
+
+ 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);
+
+ std::map<std::size_t, std::size_t> nodeIndices;
+ for(std::size_t k = 0; k < nodeTags.size(); k++)
+ nodeIndices[nodeTags[k]] = k;
+
+ #pragma omp parallel shared(nbDirichletElements, elementTags, nodeTags, nodalDisplacements, nodeIndices, nodeCoords, bcType, bcValue)
+ {
+ std::vector<elementStruct> localDirichletElementVector;
+ std::vector<elementStruct> localNeumannElementVector;
+
+ #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;
+
+ 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);
+
+ if(bcType)
+ {
+ element.bcValue[0] = bcValue;
+ localDirichletElementVector.push_back(element);
+ }
+ else
+ {
+ element.bcValue[1] = bcValue;
+ localNeumannElementVector.push_back(element);
+ }
+ }
+
+ #pragma omp critical
+ {
+ elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(), localDirichletElementVector.end());
+ elementVector.insert(elementVector.end(), localNeumannElementVector.begin(), localNeumannElementVector.end());
+ nbDirichletElements += localDirichletElementVector.size();
+ }
+ }
+ }
+ }
+ }
+
+ return nbDirichletElements;
+}
+
+
+void computeGeometricParam(elementStruct &element,
+ std::map<std::size_t, std::size_t> &nodeIndices,
+ const std::vector<double> &nodeCoords,
+ std::pair<double, double> nodalDisp1,
+ 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);
+}
+
+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;
+ #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).
+ 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];
+ }
+ }
+ }
+}
+
+void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
+ std::vector<elementStruct> &elementVector,
+ const Eigen::MatrixXd &x,
+ const std::size_t nbDirichletElements,
+ const std::string domainName)
+{
+ // Retrieving the dielectric permittivity of the BEM domain
+ std::vector<std::string> keys;
+ gmsh::onelab::getNames(keys, "(Materials).+");
+
+ double epsilon = 8.8541878128e-12; // dielectric permittivity
+ bool setEpsilon = false;
+ 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)
+ {
+ if(words[2].compare("Epsilon") == 0 && words[1].compare(domainName) == 0)
+ {
+ epsilon = value[0];
+ setEpsilon = true;
+ }
+ }
+ }
+ if(!setEpsilon)
+ std::cout << "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";
+
+ // Calculates the electrostatic pressure per element in order to provide these results to the FEM code that will handle the rest of the algorithm
+ #pragma omp parallel for schedule(guided)
+ for(size_t i = 0; i < elementVector.size(); i++)
+ {
+ if(i < nbDirichletElements)
+ elementVector[i].bcValue[1] = x(i);
+ else
+ elementVector[i].bcValue[0] = x(i);
+
+ 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;
+ }
+}
+
+
+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)
+{
+ // 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;
+}
\ No newline at end of file
diff --git a/srcs/BEM/computations.cpp b/srcs/BEM/computations.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2880c61e4731ebbcf5b6a36899893971cdb34a91
--- /dev/null
+++ b/srcs/BEM/computations.cpp
@@ -0,0 +1,334 @@
+#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>
+
+
+/*
+ * Computations
+ */
+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.hpp b/srcs/BEM/functionsBEM.hpp
index e83bb7925d3012b419fda8da266def9cb7f290b0..f8a85fbe9de02a77097217ce9bb5c251422a6d72 100644
--- a/srcs/BEM/functionsBEM.hpp
+++ b/srcs/BEM/functionsBEM.hpp
@@ -11,41 +11,167 @@
// 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, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int & nbViews, bool postProcessing, const int iteration);
-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);
+// 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).
+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
+ */
+void getBC(std::map<std::string, std::pair<bool, double>> &bcMap,
+ const std::string domainName);
+
+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);
+
+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 computeGeometricParam(elementStruct &element,
+ std::map<std::size_t, std::size_t> &nodeIndices,
+ const std::vector<double> &nodeCoords,
+ std::pair<double, double> nodalDisp1,
+ std::pair<double, double> nodalDisp2);
+
+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);
+
+void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
+ std::vector<elementStruct> &elementVector,
+ const Eigen::MatrixXd &x,
+ const std::size_t nbDirichletElements,
+ const std::string domainName);
+
+// 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);
+
+// Post processing
+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<size_t> &domainBoundaryElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+ const std::vector<elementStruct> &elementVector,
+ const int domainPhysicalGroupTag);
+
+void dataInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &data,
+ 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>>>> &data,
+ const std::vector<std::vector<std::vector<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>>>> &data,
+ 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 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);
-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);
+void displayTime(const std::chrono::time_point<std::chrono::system_clock> start, const std::chrono::time_point<std::chrono::system_clock> end, std::string function_name, 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);
\ No newline at end of file
diff --git a/srcs/BEM/mainBEM.cpp b/srcs/BEM/mainBEM.cpp
index ed6feee551bee5fb6bd5c8b5752a3fb580deb08f..e149de2fa332063a0b2ec4515b541c82ad1693a9 100644
--- a/srcs/BEM/mainBEM.cpp
+++ b/srcs/BEM/mainBEM.cpp
@@ -21,20 +21,10 @@ int main(int argc, char **argv)
int nbViews = 0;
- std::map<int,std::pair<double,double>> boundaryDisplacementMap;
+ std::map<int,std::pair<double,double>> nodalDisplacements;
//solverBEM(argc, argv, finalElementTags, electrostaticPressure, nbViews);
- solverBEM(electrostaticPressure, boundaryDisplacementMap, nbViews, true, 1); // 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, true, 1, false); // iteration number randomly set to 1
gmsh::fltk::run();
diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eb5c3b73c2e87c8658b585a6abb5609aff67b7d7
--- /dev/null
+++ b/srcs/BEM/postProcessing.cpp
@@ -0,0 +1,413 @@
+#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>
+
+
+
+/*
+ * Post processing
+ */
+
+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<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>>>> &data,
+ 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;
+ dataInitialization(data, 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 = data[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)
+ {
+ data[entityTagIndex][elementTypeIndex][elementTagIndex][j] = elementVector[commonBoundaryElementIndices[0]].bcValue[0];
+ }
+ else if(commonBoundaryElementIndices.size() == 2)
+ {
+ double phi = 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++;
+ }
+ phi += 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)
+ {
+ phi = elementVector[adjacentElementIndex].bcValue[0];
+ }
+ else
+ {
+ phi += elementVector[adjacentElementIndex].bcValue[0];
+ phi /= 2;
+ }
+ }
+ else
+ {
+ if(count == 0)
+ {
+ phi += elementVector[adjacentElementIndex].bcValue[0];
+ phi /= 2;
+ }
+ }
+
+ data[entityTagIndex][elementTypeIndex][elementTagIndex][j] = phi;
+ }
+ }
+ }
+ }
+
+ //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(data, domainElementTags, domainElementNodeTags, internalNodalPhi);
+}
+
+
+void dataInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &data,
+ 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)
+{
+ data.resize(domainElementTypes.size());
+ for(std::size_t i = 0; i < domainElementTypes.size(); i++)
+ {
+ data[i].resize(domainElementTypes[i].size());
+
+ #pragma omp parallel for schedule(guided) shared(data, 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);
+
+ data[i][j].resize(domainElementTags[i][j].size());
+ for(std::size_t k = 0; k < domainElementTags[i][j].size(); k++)
+ {
+ data[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<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>>>> &data,
+ const std::vector<std::vector<std::vector<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(data[i][j][k][l] == HUGE_VAL && internalNodalPhi.count(elementNodeTag))
+ data[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 062645f4cd3cc659c809c30b7860b8c9f0f6c0a9..f2f9c40a1368ada3f3996cae275dfaccc3e49136 100644
--- a/srcs/BEM/solverBEM.cpp
+++ b/srcs/BEM/solverBEM.cpp
@@ -11,13 +11,53 @@
#include <omp.h>
#include <list>
-void solverBEM(std::map<int, double> &electrostaticPressure, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int & nbViews, bool postProcessing, const int iteration)
+// 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).
+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)
{
+ 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);
+ }
+
+
gmsh::vectorpair physicalGroupsDimTags;
gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
- int nbDomains = 0;
+ std::vector<std::string> domainNames;
for(std::size_t i = 0; i < physicalGroupsDimTags.size(); i++)
{
@@ -32,326 +72,284 @@ void solverBEM(std::map<int, double> &electrostaticPressure, std::map<int,std::p
while(getline(ss, word, '_'))
words.push_back(word);
- if(words.size() == 3 && words[0].compare("BEM") == 0)
- {
- if(words[1].compare("domain") == 0)
- {
- if(stoi(words[2]) > nbDomains)
- {
- nbDomains = stoi(words[2]);
- }
- }
- }
+ if(words.size() == 3 && words[0].compare("BEM") == 0 && words[1].compare("domain") == 0)
+ domainNames.push_back(words[2]);
}
- if (nbDomains == 0)
+ if (domainNames.size() == 0)
{
std::cerr << "Wrong names for the BEM domains. Usage: <BEM_domain_'i'>\n";
exit(0);
}
- for(int i = 1; i <= nbDomains; i++)
- {
- singleDomainSolverBEM(electrostaticPressure, boundaryDisplacementMap, nbViews, postProcessing, iteration, i);
- }
-}
-
-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 = 1; //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.
-
- // 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);
+ int phiViewTag = gmsh::view::add("Phi");
+ int electricFieldViewTag = gmsh::view::add("Electric Field");
+ bool computePhi = true;
+ bool computeElectricField = true;
- for(size_t i = 0; i < allNodeTags.size(); i++)
+ if(iteration == 1)
{
- allNodeIndex[allNodeTags[i]] = i;
+ std::cout << "The model contains " << domainNames.size() << " BEM domain(s) \n";
+ if(postProcessing)
+ {
+ if(computePhi)
+ nbViews++;
+ if(computeElectricField)
+ nbViews++;
+ }
}
- // Print the coordinates of all nodes if dispNodeCoord != 0
- if(dispNodeCoord)
- {
- dispNodeCoordFun(allNodeTags, allCoord);
- }
+ for(std::size_t i = 0; i < domainNames.size(); i++)
+ singleDomainSolverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, iteration, "BEM_domain_" + domainNames[i], phiViewTag, electricFieldViewTag, computePhi, computeElectricField);
+}
- // Initially, these are 0 and are incremented each time a dirichlet/neumann-element is encountered.
- int dirichletCount = 0;
+// 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);
+// }
+
+// gmsh::vectorpair physicalGroupsDimTags;
+// gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
+
+// int nbDomains = 0;
+
+// for(std::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);
+
+// std::stringstream ss(name);
+// std::string word;
+// std::vector<std::string> words;
+// while(getline(ss, word, '_'))
+// words.push_back(word);
+
+// if(words.size() == 3 && words[0].compare("BEM") == 0)
+// {
+// if(words[1].compare("domain") == 0)
+// {
+// if(stoi(words[2]) > nbDomains)
+// {
+// nbDomains = stoi(words[2]);
+// }
+// }
+// }
+// }
+
+// if(iteration == 1)
+// {
+// std::cout << "The model contains " << nbDomains << " BEM domain(s) \n";
+// }
+
+// if (nbDomains == 0)
+// {
+// std::cerr << "Wrong names for the BEM domains. Usage: <BEM_domain_'i'>\n";
+// exit(0);
+// }
+
+// //#pragma omp parallel for
+// for(int i = 1; i <= nbDomains; i++)
+// {
+// singleDomainSolverBEM(electrostaticPressure, boundaryDisplacementMap, nbViews, postProcessing, iteration, i);
+// }
+// }
+
+
+
+
+
+
+// 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.
- //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);
+ 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 = false; // 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.
- //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;
+ // Used for displaying the time of different sections.
+ std::chrono::time_point<std::chrono::system_clock> start, end;
+ std::chrono::time_point<std::chrono::system_clock> startTotal, endTotal;
+ startTotal = std::chrono::system_clock::now();
- // model::mesh::getElements(elementTypes2d,elementTags2d,nodeTags2d,2);
+ // 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'.
+ // Fills the vector of boundary elements.
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);
- 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;
- }
- }
+ std::map<std::string, std::pair<int, int>> physicalGroups;
+ std::size_t nbDirichletElements = fillElementVector(elementVector, physicalGroups, nodalDisplacements, domainName);
- std::vector<std::string> names;
- gmsh::model::list(names);
- //std::cout << "the file contains " << names.size() << " model names\n";
+ end = std::chrono::system_clock::now();
+ displayTime(start, end, "'fillElementVector'", showTime);
- // Declaration of matrices A, B and vector c (the matrix type is perhaps not the most efficient).
+ // Constructs the linear system
+ start = std::chrono::system_clock::now();
+
+ // 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;
-
- 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]);
- }
- }
-
- // 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);
+ displayTime(start, end, "'fillSystem'", showTime);
-
- //If we want to display the time it takes for eigen to solve the linear system
+ // Solves 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);
- // Matrices verification
- if(dispMatrices)
- {
- dispMatricesFun(A, B, c, b, x);
- }
+ end = std::chrono::system_clock::now();
+ 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 = std::chrono::system_clock::now();
+ computeElectrostaticPressure(electrostaticPressure, elementVector, x, nbDirichletElements, domainName);
+ end = std::chrono::system_clock::now();
+ displayTime(start, end, "'computeElectrostaticPressure'", showTime);
- if(dispValues)
+ // Visualization
+ if(postProcessing)
{
- dispValuesFun(elementVector);
- }
+ // Retrieves everything related to domain of interrest.
+ start = std::chrono::system_clock::now();
- /*
- * 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;
+ 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;
- //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);
+ int domainPhysicalGroupTag = getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainEntityIndices, domainElementTypes, domainElementTags, domainElementNodeTags, physicalGroups, domainName);
- //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);
+ displayTime(start, end, "'getModel'", showTime);
+
+ // The first index indicates the "elementType", the second the "elementTag" and the third, 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;
- int viewNodeTag = gmsh::view::add("phi");
+ std::vector<std::string> names;
+ gmsh::model::list(names);
- for(size_t i = 0; i < data.size(); i++)
+ // Visualisation by ElementNodeData
+ if(visualisation)
{
- gmsh::view::addModelData(viewNodeTag,0,names[0],"ElementNodeData",elementTags2d_domain_BEM[i],data[i],0.0,1);
- }
- gmsh::view::write(viewNodeTag, "results.msh");
+ // 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);
- 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);
- nbViews++; // number of views generated (useful for FEM display of results)
+ // Computes the value of the potential at each node strictly included in the domain and stores it in the "internalNodalPhi" map
+ start = std::chrono::system_clock::now();
- 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);
+ std::map<int, double> internalNodalPhi;
+ computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector, integrationType, boundaryNodeToElements, solutionType);
+
+ end = std::chrono::system_clock::now();
+ displayTime(start, end, "'computeInternalPhi'", showTime);
+ // Associates the value of the potential to each node of an element.
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);
+ elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector, internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements);
end = std::chrono::system_clock::now();
- function_name = "compute Electric Field";
- display_time(start, end, function_name, show_time);
- 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);
- 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);
-
- if(electric_field)
- {
+ displayTime(start, end, "'elementNodeData'", showTime);
+
+ // Computes the electric field (with "ElementData" type).
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);
+ computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, false, computeElectricField);
end = std::chrono::system_clock::now();
- function_name = "compute Electric Field";
- display_time(start, end, function_name, show_time);
- 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);
- nbViews++;
- }
+ displayTime(start, end, "'computeElementData' (only the electric field)", showTime);
- /*
- * Convergence
- */
-
- if(convergence)
+ // Displays what is needed.
+ if(computePhi)
+ displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi, nbViews);
+ if(computeElectricField)
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, nbViews+1);
+ }
+ // Visualisation by ElementData
+ else
{
- convergenceFun(entities_tags_2d_domain_BEM,tags,phi);
+ // Computes the potential and the electric field.
+ start = std::chrono::system_clock::now();
+ computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, computePhi, computeElectricField);
+ end = std::chrono::system_clock::now();
+ displayTime(start, end, "'computeElementData'", showTime);
+
+ // Displays what is needed.
+ if(computePhi)
+ displayData(phiViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalPhi, nbViews);
+ if(computeElectricField)
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, nbViews+1);
+
+ // Study of the error on a precise model.
+ if(convergence)
+ {
+ std::string convergenceModelName = "BEM_convergence";
+ if(names[0].compare(convergenceModelName) == 0)
+ {
+ std::cout << "\nNumber of boundary elements: " << elementVector.size() << "\n";
+ convergenceFun(domainEntityTags, domainElementTypes, domainElementTags, elementalPhi, elementalElectricField);
+ }
+ else
+ std::cout << "\n\nWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '" << convergenceModelName << "' model if you want to check for convergence.\n\n";
+ }
}
-
- //For ease of use, a function handles the code for the display.
- displayData("phi", names[0], "results.msh", tags, phi);
- nbViews++; // number of views generated (useful for FEM display of results)
- }
-
- //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;
- for(size_t i = 0; i < elementVector.size(); i++)
- {
- 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 time of the code BEM";
- display_time(start_total, end_total, function_name, show_time);
+ endTotal = std::chrono::system_clock::now();
+ 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..68528f03d7ff644487b48d4fe952f7b1da430933
--- /dev/null
+++ b/srcs/BEM/verification.cpp
@@ -0,0 +1,117 @@
+#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>
+
+/*
+ * Verification
+ */
+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 displayTime(const std::chrono::time_point<std::chrono::system_clock> start,
+ const std::chrono::time_point<std::chrono::system_clock> end,
+ std::string function_name, bool &showTime)
+{
+ //displays the time of the desired code section
+ if(showTime)
+ {
+ 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 << "Relative error on the electric potential: " << errorPhi << "\n";
+ std::cout << "Relative error on the electric field: " << errorE << "\n";
+ std::cout << "Relative error on the electric field (without the edges of the domain): " << errorEWindow << "\n\n";
+}
\ No newline at end of file