diff --git a/srcs/BEM/boundaryElements.cpp b/srcs/BEM/boundaryElements.cpp
index 7884bf59fbe7fa6c57c1b97b7fb2db77eff90e99..bdfa2d7408deecb5ec785e2be9e745c6d5fc9eb2 100644
--- a/srcs/BEM/boundaryElements.cpp
+++ b/srcs/BEM/boundaryElements.cpp
@@ -12,12 +12,12 @@
#include <algorithm>
#include <list>
-void getBC(std::map<std::string, std::pair<bool, double>> &bcMap,
+// 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)
{
- // 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
+ // Extracts all the boundary conditions from "ONELAB".
std::vector<std::string> bcKeys;
gmsh::onelab::getNames(bcKeys, "Boundary Conditions.+");
@@ -32,26 +32,31 @@ void getBC(std::map<std::string, std::pair<bool, double>> &bcMap,
while(getline(ss, word, '/'))
words.push_back(word);
- // Only focus on the BEM boundary conditions.
+ // Only focus on the right BEM domain.
if(words.size() == 4 && words[2].compare(domainName) == 0)
{
if(words[3].compare("dirichlet") == 0)
- bcMap[words[1]] = {true, values[0]};
+ boundaryConditions[words[1]] = {true, values[0]};
else if(words[3].compare("neumann") == 0)
- bcMap[words[1]] = {false, values[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::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::size_t nbDirichletElements = 0; // Will be incremented later
- std::map<std::string, std::pair<bool, double>> bcMap;
- getBC(bcMap, domainName);
+ 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;
@@ -69,43 +74,49 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
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())
+ if(boundaryConditions.find(physicalGroupName) == boundaryConditions.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;
+ 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};
+ 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++)
{
@@ -119,7 +130,7 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
if(nbNodesPerElement == 3)
element.midNodeTag = nodeTags[nbNodesPerElement*el + 2];
else
- element.midNodeTag = 0;
+ 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]))
@@ -128,6 +139,7 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
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;
@@ -140,6 +152,8 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
}
}
+ // 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());
@@ -154,12 +168,15 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
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,
- std::pair<double, double> nodalDisp1,
- std::pair<double, double> nodalDisp2)
+ 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];
@@ -180,6 +197,10 @@ void computeGeometricParam(elementStruct &element,
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,
@@ -190,6 +211,8 @@ void fillSystem(Eigen::MatrixXd &A,
{
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++)
{
@@ -213,7 +236,7 @@ void fillSystem(Eigen::MatrixXd &A,
g = computeAnalyticalG(elementVector[i].midX, elementVector[i].midY, elementVector[j]);
}
- // Storage of h and g at the right place (in matrices A and B).
+ // Storage of h and g at the right place (in matrices A and B and vector c).
if(j < nbDirichletElements)
{
A(i,j) = g;
@@ -230,13 +253,18 @@ void fillSystem(Eigen::MatrixXd &A,
}
}
+// 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)
{
- // Retrieving the dielectric permittivity of the BEM domain
+ // 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).+");
@@ -244,11 +272,10 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
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
+
+ // Expected key structure is "type/group_name/field"
std::stringstream ss(key);
std::vector<std::string> words;
std::string word;
@@ -264,83 +291,21 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
}
}
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";
+ 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";
- // 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++)
{
+ // 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;
}
-}
-
-
-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
index 2880c61e4731ebbcf5b6a36899893971cdb34a91..c35b66d0c76ecebd80aff3deae8237ecc8571920 100644
--- a/srcs/BEM/computations.cpp
+++ b/srcs/BEM/computations.cpp
@@ -12,10 +12,6 @@
#include <algorithm>
#include <list>
-
-/*
- * Computations
- */
double computeAnalyticalH(const double x,
const double y,
const elementStruct &element)
diff --git a/srcs/BEM/functionsBEM.hpp b/srcs/BEM/functionsBEM.hpp
index f8a85fbe9de02a77097217ce9bb5c251422a6d72..e705a73356775bcb6d84184ecb05b2d99d7b7f7a 100644
--- a/srcs/BEM/functionsBEM.hpp
+++ b/srcs/BEM/functionsBEM.hpp
@@ -30,7 +30,8 @@ struct elementStruct{
// 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).
+// iteration). "untangleMesh" unstangles every BEM domain if set to 'true' (in case the "nodalDisplacements" map contains
+// non-zero displacements).
void solverBEM(std::map<int, double> &electrostaticPressure,
int & nbViews,
std::map<int,std::pair<double,double>> & nodalDisplacements,
@@ -57,31 +58,36 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
/*
* Boundary element
*/
-void getBC(std::map<std::string, std::pair<bool, double>> &bcMap,
+
+// Fills the "boundaryConditions" map associating to each physical group (a string) the type of boundary condition ('true':
+// dirichlet, 'false': neumann) and its value (if any). "domainName" is the name of the current BEM domain.
+void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
const std::string domainName);
+// Fills the "elementVector" vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the
+// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also fills the
+// "physicalGroups" map associating to each (physical group) name, its tag and dimension. "domainName" is the name of the current
+// BEM domain. "nodalDisplacements" is a map associating the (x,y) displacement to a serie of node tags. Returns the number of
+// 'dirichlet elements' (useful for the construction/resolution of the system).
std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
- std::map<std::string, std::pair<int, int>> &physicalGroups,
- std::map<int, std::pair<double, double>> &nodalDisplacements,
- const std::string domainName);
-
-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);
-
+ std::map<std::string, std::pair<int, int>> &physicalGroups,
+ std::map<int, std::pair<double, double>> &nodalDisplacements,
+ const std::string domainName);
+
+// Computes the useful x of a given element and stores it in "element". "nodeIndinces" is a map associating the index of any
+// node tag for accessing its coordinates in the "nodeCoords" vector. "nodalDisp1" and "nodalDisp2" are pairs representing the
+// potential displacement of the first and second node (respectivelly) of the element. If an element has 3 nodes, the first and
+// second nodes are considered to be the extremities of that element.
void computeGeometricParam(elementStruct &element,
std::map<std::size_t, std::size_t> &nodeIndices,
const std::vector<double> &nodeCoords,
- std::pair<double, double> nodalDisp1,
- std::pair<double, double> nodalDisp2);
+ const std::pair<double, double> &nodalDisp1,
+ const std::pair<double, double> &nodalDisp2);
+// Fills the "A", "B", "c" matrices (/vector) that represent the linear system. "elementVector" is a vector storing the boundary
+// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "nbDirichletElements" is the number of
+// 'dirichlet elements'. "solutionType" indicates whether the elements of the matrices are computed analytically ('true') or
+// numerically ('false'). "integrationType" is the type of numerical integration used (in case of numerical computation).
void fillSystem(Eigen::MatrixXd &A,
Eigen::MatrixXd &B,
Eigen::MatrixXd &c,
@@ -90,26 +96,31 @@ void fillSystem(Eigen::MatrixXd &A,
const std::string integrationType,
const bool solutionType);
+// Fills the "electrostaticPressure" map associating each (boundary) element tag with its electrostatic pressure. Also sets the
+// unknown of each element (either the potential of the electric field) of the "elementVector" vector storing the boundary
+// elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM domain. "x" is the solution of the linear
+// system. "nbDirichletElements" is the number of 'dirichlet elements', "domainName" is the name of the considered BEM domain.
void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
std::vector<elementStruct> &elementVector,
const Eigen::MatrixXd &x,
const std::size_t nbDirichletElements,
const std::string domainName);
-// 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
+ */
+
+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);
-// Post processing
void computeInternalPhi(std::map<int, double> &internalNodalPhi,
const std::vector<double> &domainNodeCoords,
const std::vector<size_t> &domainNodeTags,
@@ -122,23 +133,23 @@ void computeInternalPhi(std::map<int, double> &internalNodalPhi,
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,
+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<std::vector<std::vector<std::size_t>>> &domainElementTags,
const std::vector<elementStruct> &elementVector,
const int domainPhysicalGroupTag);
-void dataInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &data,
+void phiInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
std::map<std::size_t, std::vector<std::size_t>> &position,
const std::vector<std::vector<int>> &domainElementTypes,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags);
-void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &data,
- const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
std::map<int, double> &internalNodalPhi);
-void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>> &data,
+void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
const std::vector<std::vector<int>> &domainElementTypes,
std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
@@ -166,12 +177,33 @@ void displayData(const int viewTag,
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 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);
+/*
+ * Verification
+ */
+
+void displayTime(const double start,
+ const double end,
+ std::string functionName, const bool showTime);
void convergenceFun(const std::vector<int> &domainEntityTags,
const std::vector<std::vector<int>> &domainElementTypes,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalPhi,
- const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalElectricField);
\ No newline at end of file
+ const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalElectricField);
+
+
+/*
+ * Computations
+ */
+
+double computeAnalyticalH(const double x, const double y, const elementStruct &element);
+double computeNumericalG(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradHx(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradHy(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradGx(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeNumericalGradGy(const double x, const double y, const elementStruct &element, const std::string integrationType);
+double computeAnalyticalGradHx(const double x, const double y, const elementStruct &element);
+double computeAnalyticalGradHy(const double x, const double y, const elementStruct &element);
+double computeAnalyticalGradGx(const double x, const double y, const elementStruct &element);
+double computeAnalyticalGradGy(const double x, const double y, const elementStruct &element);
+double computeAnalyticalG(const double x, const double y, const elementStruct &element);
\ No newline at end of file
diff --git a/srcs/BEM/mainBEM.cpp b/srcs/BEM/mainBEM.cpp
index e149de2fa332063a0b2ec4515b541c82ad1693a9..35655adb9b513d22220df01b905da9a8178b7b90 100644
--- a/srcs/BEM/mainBEM.cpp
+++ b/srcs/BEM/mainBEM.cpp
@@ -18,15 +18,15 @@ int main(int argc, char **argv)
// will map the 1d element tag onto the corresponding electroPressure resulting from BEM code
std::map<int, double> electrostaticPressure;
-
int nbViews = 0;
-
std::map<int,std::pair<double,double>> nodalDisplacements;
+ const bool postProcessing = true;
+ const bool meshUntangler = false;
+
+ solverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, 1, meshUntangler); // iteration number randomly set to 1
- //solverBEM(argc, argv, finalElementTags, electrostaticPressure, nbViews);
- solverBEM(electrostaticPressure, nbViews, nodalDisplacements, true, 1, false); // iteration number randomly set to 1
-
- gmsh::fltk::run();
+ if(postProcessing)
+ gmsh::fltk::run();
gmsh::finalize();
return 0;
diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp
index eb5c3b73c2e87c8658b585a6abb5609aff67b7d7..c5854559391eff43301a958a05c5861d9c752c03 100644
--- a/srcs/BEM/postProcessing.cpp
+++ b/srcs/BEM/postProcessing.cpp
@@ -12,11 +12,68 @@
#include <algorithm>
#include <list>
+int getModel(std::vector<std::size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ std::vector<double> &domainNodeCoords,
+ std::vector<int> &domainEntityTags,
+ std::map<int, std::size_t> &domainEntityIndices,
+ std::vector<std::vector<int>> &domainElementTypes,
+ std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
+ std::map<std::string, std::pair<int, int>> &physicalGroups,
+ const std::string domainName)
+{
+ // Retrieves everything only for the domain of interrest.
+ int physicalGroupDim = physicalGroups[domainName].first;
+ int physicalGroupTag = physicalGroups[domainName].second;
+
+ gmsh::model::mesh::getNodesForPhysicalGroup(physicalGroupDim, physicalGroupTag, domainNodeTags, domainNodeCoords);
+
+ // The "domainNodeIndices" map is used for finding the position of the coordinates of a node in "domainNodeCoords"
+ // from its tag.
+ for(std::size_t index = 0; index < domainNodeTags.size(); index++)
+ domainNodeIndices[domainNodeTags[index]] = index;
+
+ gmsh::model::getEntitiesForPhysicalGroup(physicalGroupDim, physicalGroupTag, domainEntityTags);
+ domainElementTypes.resize(domainEntityTags.size());
+ domainElementTags.resize(domainEntityTags.size());
+ for(std::size_t j = 0; j < domainEntityTags.size(); j++)
+ {
+ domainEntityIndices[domainEntityTags[j]] = j;
+
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> nodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, physicalGroupDim, domainEntityTags[j]);
+
+ domainElementTypes[j].resize(elementTypes.size());
+ domainElementTags[j].resize(elementTypes.size());
+ for(std::size_t k = 0; k < elementTypes.size(); k++)
+ {
+ int nbNodes = nodeTags[k].size() / elementTags[k].size();
+ std::vector<std::size_t> elementNodeTags(nbNodes);
+
+ domainElementTypes[j][k] = elementTypes[k];
+ domainElementTags[j][k].resize(elementTags[k].size());
+
+ for(std::size_t l = 0; l < elementTags[k].size(); l++)
+ {
+ domainElementTags[j][k][l] = elementTags[k][l];
+ int n = 0;
+ for(std::size_t m = l*nbNodes; m < (l+1)*nbNodes; m++)
+ {
+ elementNodeTags[n] = nodeTags[k][m];
+ n++;
+ }
+
+ domainElementNodeTags[elementTags[k][l]] = elementNodeTags;
+ }
+ }
+ }
+ return physicalGroupTag;
+}
-/*
- * Post processing
- */
void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
const std::vector<elementStruct> &elementVector)
@@ -39,13 +96,13 @@ void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boun
}
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)
+ const std::vector<double> &domainNodeCoords,
+ const std::vector<std::size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ const std::vector<elementStruct> &elementVector,
+ const std::string integrationType,
+ std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
+ const bool &solutionType)
{
double x = 0;
double y = 0;
@@ -79,12 +136,10 @@ void computeInternalPhi(std::map<int, double> &internalNodalPhi,
}
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,
+void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
const std::vector<std::vector<int>> &domainElementTypes,
std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
@@ -96,7 +151,7 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
// 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);
+ phiInitialization(phi, position, domainElementTypes, domainElementTags);
// // The first part of the algorithm consists in recovering the tag of the 2d elements of the domain which have at least one node in common with the boundary
// // We initialise the first dimension of the vector that will contain all 2d elements in common with the boundary with the number of element types present in the domain.
@@ -112,7 +167,7 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
int entityTagIndex = position[elementTag][0];
int elementTypeIndex = position[elementTag][1];
int elementTagIndex = position[elementTag][2];
- int nbNodes = data[entityTagIndex][elementTypeIndex][elementTagIndex].size();
+ int nbNodes = phi[entityTagIndex][elementTypeIndex][elementTagIndex].size();
std::vector<std::size_t> nodeTags = domainElementNodeTags[elementTag];
// The third loop is on the nodes of a specific elementTag in 2d in common with the boundary.
@@ -127,11 +182,11 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
// 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];
+ phi[entityTagIndex][elementTypeIndex][elementTagIndex][j] = elementVector[commonBoundaryElementIndices[0]].bcValue[0];
}
else if(commonBoundaryElementIndices.size() == 2)
{
- double phi = 0.0;
+ double currentPhi = 0.0;
int count = 0;
// First element.
@@ -147,7 +202,7 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
{
count++;
}
- phi += elementVector[adjacentElementIndex].bcValue[0];
+ currentPhi += elementVector[adjacentElementIndex].bcValue[0];
// Second element.
@@ -163,45 +218,45 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
{
if(count == 0)
{
- phi = elementVector[adjacentElementIndex].bcValue[0];
+ currentPhi = elementVector[adjacentElementIndex].bcValue[0];
}
else
{
- phi += elementVector[adjacentElementIndex].bcValue[0];
- phi /= 2;
+ currentPhi += elementVector[adjacentElementIndex].bcValue[0];
+ currentPhi /= 2;
}
}
else
{
if(count == 0)
{
- phi += elementVector[adjacentElementIndex].bcValue[0];
- phi /= 2;
+ currentPhi += elementVector[adjacentElementIndex].bcValue[0];
+ currentPhi /= 2;
}
}
- data[entityTagIndex][elementTypeIndex][elementTagIndex][j] = phi;
+ phi[entityTagIndex][elementTypeIndex][elementTagIndex][j] = currentPhi;
}
}
}
}
//The last step is to loop over all the 2d elements in the domain. If the HUGE_VAL value is still present, this means that the element node does not belong to the boundary (since the value of the potential has not been modified) and is therefore part of the domain. We then use the map which associates the value of the potential with the tag of the node inside the domain.
- boundaryElementFinalization(data, domainElementTags, domainElementNodeTags, internalNodalPhi);
+ boundaryElementFinalization(phi, 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)
+void phiInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ std::map<std::size_t, std::vector<std::size_t>> &position,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<size_t>>> &domainElementTags)
{
- data.resize(domainElementTypes.size());
+ phi.resize(domainElementTypes.size());
for(std::size_t i = 0; i < domainElementTypes.size(); i++)
{
- data[i].resize(domainElementTypes[i].size());
+ phi[i].resize(domainElementTypes[i].size());
- #pragma omp parallel for schedule(guided) shared(data, domainElementTypes, domainElementTags, position)
+ #pragma omp parallel for schedule(guided) shared(phi, domainElementTypes, domainElementTags, position)
for(std::size_t j = 0; j < domainElementTypes[i].size(); j++)
{
std::string elementName;
@@ -212,10 +267,10 @@ void dataInitialization(std::vector<std::vector<std::vector<std::vector<double>>
int nbPrimaryNodes;
gmsh::model::mesh::getElementProperties(domainElementTypes[i][j], elementName, dim, order, nbNodes, localNodeCoord, nbPrimaryNodes);
- data[i][j].resize(domainElementTags[i][j].size());
+ phi[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);
+ phi[i][j][k].resize(nbNodes, HUGE_VAL);
std::vector<std::size_t> pos = {i, j, k};
#pragma omp critical
@@ -227,7 +282,7 @@ void dataInitialization(std::vector<std::vector<std::vector<std::vector<double>>
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<std::vector<std::vector<std::size_t>>> &domainElementTags,
const std::vector<elementStruct> &elementVector,
const int domainPhysicalGroupTag)
{
@@ -273,8 +328,8 @@ void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
}
-void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &data,
- const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
+void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
std::map<int, double> &internalNodalPhi)
{
@@ -289,8 +344,8 @@ void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector
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];
+ if(phi[i][j][k][l] == HUGE_VAL && internalNodalPhi.count(elementNodeTag))
+ phi[i][j][k][l] = internalNodalPhi[elementNodeTag];
}
}
}
@@ -406,7 +461,7 @@ void displayData(const int viewTag,
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);
diff --git a/srcs/BEM/solverBEM.cpp b/srcs/BEM/solverBEM.cpp
index f2f9c40a1368ada3f3996cae275dfaccc3e49136..67ab21db5341f66df856edbd53451b30359f963e 100644
--- a/srcs/BEM/solverBEM.cpp
+++ b/srcs/BEM/solverBEM.cpp
@@ -16,7 +16,8 @@
// 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).
+// iteration). "untangleMesh" unstangles every BEM domain if set to 'true' (in case the "nodalDisplacements" map contains
+// non-zero displacements).
void solverBEM(std::map<int, double> &electrostaticPressure,
int & nbViews,
std::map<int,std::pair<double,double>> & nodalDisplacements,
@@ -24,40 +25,15 @@ void solverBEM(std::map<int, double> &electrostaticPressure,
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]);
+ bool computePhi = true; // true: computes the potential.
+ bool computeElectricField = true; // true: computes the electric field.
- gmsh::model::mesh::optimize("UntangleMeshGeometry", false, 1, dimTagsBis);
- }
-
-
+ // Retrieves the names of the BEM domains as well as their physical dim/tag.
gmsh::vectorpair physicalGroupsDimTags;
gmsh::model::getPhysicalGroups(physicalGroupsDimTags);
std::vector<std::string> domainNames;
+ std::vector<std::pair<int, int>> domainPhysicalGroupDimTags;
for(std::size_t i = 0; i < physicalGroupsDimTags.size(); i++)
{
@@ -73,118 +49,71 @@ void solverBEM(std::map<int, double> &electrostaticPressure,
words.push_back(word);
if(words.size() == 3 && words[0].compare("BEM") == 0 && words[1].compare("domain") == 0)
+ {
domainNames.push_back(words[2]);
+ domainPhysicalGroupDimTags.push_back({physicalGroupsDim, physicalGroupsTag});
+ }
}
+ // Useful for dealing with the untangler.
+ if(untangleMesh)
+ // Untangles each BEM domain.
+ for(std::size_t i = 0; i < domainNames.size(); i++)
+ {
+ int dim = domainPhysicalGroupDimTags[i].first;
+ int tag = domainPhysicalGroupDimTags[i].second;
+
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+ gmsh::vectorpair dimTagsBis(tags.size());
+ for(size_t i = 0; i < tags.size(); i++)
+ dimTagsBis[i] = std::pair<int,int>(dim, tags[i]);
+
+ gmsh::model::mesh::optimize("UntangleMeshGeometry", false, 1, dimTagsBis);
+
+ }
+
+
if (domainNames.size() == 0)
{
- std::cerr << "Wrong names for the BEM domains. Usage: <BEM_domain_'i'>\n";
+ std::cerr << "Error: No BEM domain found ! Usage: <BEM_domain_'domainName'>\n";
exit(0);
}
- int phiViewTag = gmsh::view::add("Phi");
- int electricFieldViewTag = gmsh::view::add("Electric Field");
- bool computePhi = true;
- bool computeElectricField = true;
-
if(iteration == 1)
{
- std::cout << "The model contains " << domainNames.size() << " BEM domain(s) \n";
- if(postProcessing)
- {
- if(computePhi)
- nbViews++;
- if(computeElectricField)
- nbViews++;
- }
+ std::cout << "--------------------[BEM SOLVER]--------------------\n";
+ std::cout << "The model contains " << domainNames.size() << " BEM domain(s):\n";
+ }
+
+ // Generates viewTags only if the postProcessing is enabled.
+ int phiViewTag = 0;
+ int electricFieldViewTag = 0;
+ if(postProcessing)
+ {
+ if(computePhi)
+ phiViewTag = gmsh::view::add("Phi");
+ if(computeElectricField)
+ electricFieldViewTag = gmsh::view::add("Electric Field");
}
+ // Runs the solver for each BEM domain.
for(std::size_t i = 0; i < domainNames.size(); i++)
singleDomainSolverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, iteration, "BEM_domain_" + domainNames[i], phiViewTag, electricFieldViewTag, computePhi, computeElectricField);
-}
-
-// 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);
-// }
-// }
-
-
+ if(iteration == 1)
+ std::cout << "--------------------[BEM SOLVER]--------------------\n\n";
+ // Updates the number of views.
+ if(postProcessing)
+ {
+ if(computePhi)
+ nbViews++;
+ if(computeElectricField)
+ nbViews++;
+ }
+}
// Fills the "electrostaticPressure" map. The "nodalDisplacements", "postProcessing" and "iteration" arguments are described in
@@ -207,29 +136,33 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
bool showTime = false; // true: displays the different time measurements.
bool convergence = false; // true: computes the error between the analytical solution and the numerical one on a simple configuration (/!\ use this with the BEM_convergence.geo model).
- bool visualisation = false; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of the potential (Remark: electric field is always displayed with 'ElementData' type).
+ bool visualisation = true; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of the potential (Remark: electric field is always displayed with 'ElementData' type).
bool solutionType = true; // true: analytical computations, false: numerical computations.
// Used for displaying the time of different sections.
- 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();
+ // double start, end;
+ double startTotal, endTotal;
+ double start, end;
+ startTotal = omp_get_wtime();
// Initializes a vector of elements (that will contain all the elements of interrest)
std::vector<elementStruct> elementVector;
std::string integrationType = "Gauss10"; // Used for the numerical integrations when "solutionType" = 'false'.
// Fills the vector of boundary elements.
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
std::map<std::string, std::pair<int, int>> physicalGroups;
std::size_t nbDirichletElements = fillElementVector(elementVector, physicalGroups, nodalDisplacements, domainName);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'fillElementVector'", showTime);
+ if(iteration == 1)
+ std::cout << "\n- BEM domain '" << domainName << "': contains " << elementVector.size() << " boundary elements.\n";
+
// Constructs the linear system
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
// The matrix type is perhaps not the most efficient.
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(elementVector.size(), elementVector.size());
@@ -237,29 +170,31 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
Eigen::MatrixXd c = Eigen::MatrixXd::Zero(elementVector.size(), 1);
fillSystem(A, B, c, elementVector, nbDirichletElements, integrationType, solutionType);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'fillSystem'", showTime);
// Solves the linear system.
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
Eigen::MatrixXd b = B*c;
- Eigen::MatrixXd x = A.fullPivLu().solve(b);
+ Eigen::MatrixXd x = A.partialPivLu().solve(b);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "Solve linear system", showTime);
// Computes the electrostatic pressure on each boundary element.
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
computeElectrostaticPressure(electrostaticPressure, elementVector, x, nbDirichletElements, domainName);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'computeElectrostaticPressure'", showTime);
// Visualization
if(postProcessing)
{
+ int viewIndex = nbViews;
+
// Retrieves everything related to domain of interrest.
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
std::vector<std::size_t> domainNodeTags;
std::map<std::size_t, std::size_t> domainNodeIndices;
@@ -273,12 +208,12 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
int domainPhysicalGroupTag = getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainEntityIndices, domainElementTypes, domainElementTags, domainElementNodeTags, physicalGroups, domainName);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
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).
+ // The first index indicates the "entityTag", the second the "elementType", the third the "elementTag", and the fourth,
+ // either the nodesTags associated with the element (for elementNodalPhi) or the computed values (for "elemental...") (1
+ // for the potential and 3 for the electric field).
std::vector<std::vector<std::vector<std::vector<double>>>> elementalPhi;
std::vector<std::vector<std::vector<std::vector<double>>>> elementNodalPhi;
std::vector<std::vector<std::vector<std::vector<double>>>> elementalElectricField;
@@ -294,62 +229,66 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
computeNodeToElements(boundaryNodeToElements, elementVector);
// Computes the value of the potential at each node strictly included in the domain and stores it in the "internalNodalPhi" map
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
std::map<int, double> internalNodalPhi;
computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector, integrationType, boundaryNodeToElements, solutionType);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'computeInternalPhi'", showTime);
// Associates the value of the potential to each node of an element.
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector, internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'elementNodeData'", showTime);
// Computes the electric field (with "ElementData" type).
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, false, computeElectricField);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'computeElementData' (only the electric field)", showTime);
// Displays what is needed.
if(computePhi)
- displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi, nbViews);
+ {
+ displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi, viewIndex);
+ viewIndex++;
+ }
if(computeElectricField)
- displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, nbViews+1);
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex);
}
// Visualisation by ElementData
else
{
// Computes the potential and the electric field.
- start = std::chrono::system_clock::now();
+ start = omp_get_wtime();
computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, computePhi, computeElectricField);
- end = std::chrono::system_clock::now();
+ end = omp_get_wtime();
displayTime(start, end, "'computeElementData'", showTime);
// Displays what is needed.
if(computePhi)
- displayData(phiViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalPhi, nbViews);
+ {
+ displayData(phiViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalPhi, viewIndex);
+ viewIndex++;
+ }
if(computeElectricField)
- displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, nbViews+1);
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex);
// Study of the error on a precise model.
if(convergence)
{
+ std::cout << "\tCONVERGENCE:\n";
std::string convergenceModelName = "BEM_convergence";
if(names[0].compare(convergenceModelName) == 0)
- {
- 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";
+ std::cout << "\tWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '" << convergenceModelName + ".geo" << "' model if you want to check for convergence.\n\n";
}
}
}
- endTotal = std::chrono::system_clock::now();
+ endTotal = omp_get_wtime();
displayTime(startTotal, endTotal, "total time of the code BEM", showTime);
}
\ No newline at end of file
diff --git a/srcs/BEM/verification.cpp b/srcs/BEM/verification.cpp
index 68528f03d7ff644487b48d4fe952f7b1da430933..48d7bbde1138e9dd2ccb7bd97b1502ea3555a34e 100644
--- a/srcs/BEM/verification.cpp
+++ b/srcs/BEM/verification.cpp
@@ -12,45 +12,36 @@
#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)
+void displayTime(const double start,
+ const double end,
+ std::string functionName, const 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;
+ double duration = end - start;
+ std::cout << functionName << "\n";
+ std::cout << "Elapsed time: " << 1000000*duration << " [µs]\tor " << duration << " [s]\n\n";
+
+ // std::cout << std::endl;
+ // std::cout << function_name << std::endl;
+ // std::cout << "Elapsed time in nanoseconds: "
+ // << std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count()
+ // << " ns" << std::endl;
+
+ // std::cout << "Elapsed time in microseconds: "
+ // << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
+ // << " µs" << std::endl;
+
+ // std::cout << "Elapsed time in milliseconds: "
+ // << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
+ // << " ms" << std::endl;
+
+ // std::cout << "Elapsed time in seconds: "
+ // << std::chrono::duration_cast<std::chrono::seconds>(end - start).count()
+ // << " sec";
+ // std::cout << std::endl;
+ // std::cout << std::endl;
}
}
@@ -111,7 +102,8 @@ void convergenceFun(const std::vector<int> &domainEntityTags,
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";
+
+ std::cout << "\tRelative error on the electric potential: " << errorPhi << ".\n";
+ std::cout << "\tRelative error on the electric field: " << errorE << ".\n";
+ std::cout << "\tRelative error on the electric field (without the edges of the domain): " << errorEWindow << ".\n\n";
}
\ No newline at end of file