Py bem

srcs/BEM/boundaryElements.cpp

+#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