From 81ad07f74793487c2568c132748397835911314c Mon Sep 17 00:00:00 2001
From: Pierre-Yves <PYGoffin@student.uliege.be>
Date: Sat, 21 May 2022 02:50:52 +0200
Subject: [PATCH] comments BEM
---
srcs/BEM/boundaryElements.cpp | 47 +++--
srcs/BEM/computations.cpp | 40 +++--
srcs/BEM/functionsBEM.hpp | 267 ++++++++++++++++++++-------
srcs/BEM/mainBEM.cpp | 6 +-
srcs/BEM/postProcessing.cpp | 330 ++++++++++++++++++++--------------
srcs/BEM/solverBEM.cpp | 66 ++++---
srcs/BEM/verification.cpp | 46 ++---
7 files changed, 498 insertions(+), 304 deletions(-)
diff --git a/srcs/BEM/boundaryElements.cpp b/srcs/BEM/boundaryElements.cpp
index bdfa2d7..54dd80a 100644
--- a/srcs/BEM/boundaryElements.cpp
+++ b/srcs/BEM/boundaryElements.cpp
@@ -1,16 +1,4 @@
-#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.
@@ -44,12 +32,13 @@ void getBC(std::map<std::string, std::pair<bool, double>> &boundaryConditions,
}
// 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).
+// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also returns
+// "domainPhysicalGroupDim" and "domainPhysicalGroupTag", the dimension and the tag (respectivelly) of the BEM domain of
+// interest. "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,
+ int &domainPhysicalGroupDim,
+ int &domainPhysicalGroupTag,
std::map<int, std::pair<double, double>> &nodalDisplacements,
const std::string domainName)
{
@@ -68,11 +57,14 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
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(physicalGroupName.compare(domainName) == 0)
+ {
+ domainPhysicalGroupDim = physicalGroupDim;
+ domainPhysicalGroupTag = physicalGroupTag;
+ }
+
// If no boundary condition on a given physical curve, don't do anything.
if(boundaryConditions.find(physicalGroupName) == boundaryConditions.end())
continue;
@@ -103,7 +95,8 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
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);
+ 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;
@@ -156,8 +149,10 @@ std::size_t fillElementVector(std::vector<elementStruct> &elementVector,
// the vector.
#pragma omp critical
{
- elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(), localDirichletElementVector.end());
- elementVector.insert(elementVector.end(), localNeumannElementVector.begin(), localNeumannElementVector.end());
+ elementVector.insert(elementVector.begin(), localDirichletElementVector.begin(),
+ localDirichletElementVector.end());
+ elementVector.insert(elementVector.end(), localNeumannElementVector.begin(),
+ localNeumannElementVector.end());
nbDirichletElements += localDirichletElementVector.size();
}
}
@@ -291,7 +286,11 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
}
}
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";
+ {
+ std::cout << "\tWarning: You forgot to set the dielectric permittivity of the material for the domain '" << domainName;
+ std::cout << "' in the .geo file and it thus has been fixed to the dielectric permittivity of the vacuum: " << epsilon;
+ std::cout << " [F/m].\n";
+ }
#pragma omp parallel for schedule(guided)
for(size_t i = 0; i < elementVector.size(); i++)
diff --git a/srcs/BEM/computations.cpp b/srcs/BEM/computations.cpp
index c35b66d..9666756 100644
--- a/srcs/BEM/computations.cpp
+++ b/srcs/BEM/computations.cpp
@@ -1,17 +1,7 @@
-#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>
+// Computes the 'h' element analytically. "x" and "y" are the coordinates at which 'h' should be computed. "element" is the
+// element 'j' while will helps to compute 'h'.
double computeAnalyticalH(const double x,
const double y,
const elementStruct &element)
@@ -41,6 +31,8 @@ double computeAnalyticalH(const double x,
return alpha / (2*M_PI);
}
+// Computes the 'g' element analytically. "x" and "y" are the coordinates at which 'g' should be computed. "element" is the
+// element 'j' while will helps to compute 'g'.
double computeAnalyticalG(const double x,
const double y,
const elementStruct &element)
@@ -58,6 +50,8 @@ double computeAnalyticalG(const double x,
return g_analytique;
}
+// Computes the 'hx' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradHx' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradHx'.
double computeAnalyticalGradHx(const double x,
const double y,
const elementStruct &element)
@@ -86,6 +80,8 @@ double computeAnalyticalGradHx(const double x,
return gradHx_analytique;
}
+// Computes the 'hy' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradHy' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradHy'.
double computeAnalyticalGradHy(const double x,
const double y,
const elementStruct &element)
@@ -114,6 +110,8 @@ double computeAnalyticalGradHy(const double x,
return gradHy_analytique;
}
+// Computes the 'gx' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradGx' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradGx'.
double computeAnalyticalGradGx(const double x,
const double y,
const elementStruct &element)
@@ -133,6 +131,8 @@ double computeAnalyticalGradGx(const double x,
return gradGx_analytique;
}
+// Computes the 'gy' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradGy' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradGy'.
double computeAnalyticalGradGy(const double x,
const double y,
const elementStruct &element)
@@ -152,7 +152,9 @@ double computeAnalyticalGradGy(const double x,
return gradGy_analytique;
}
-
+// Computes the 'g' element numerically. "x" and "y" are the coordinates at which 'g' should be computed. "element" is the
+// element 'j' while will helps to compute 'g'. "integrationType" is the type of numerical integration used (in case of numerical
+// computation).
double computeNumericalG(const double x,
const double y,
const elementStruct &element,
@@ -182,6 +184,9 @@ double computeNumericalG(const double x,
return g;
}
+// Computes the 'hx' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradHx' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradHx'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
double computeNumericalGradHx(const double x,
const double y,
const elementStruct &element,
@@ -226,6 +231,9 @@ double computeNumericalGradHx(const double x,
return gradHx;
}
+// Computes the 'hy' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradHy' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradHy'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
double computeNumericalGradHy(const double x,
const double y,
const elementStruct &element,
@@ -270,6 +278,9 @@ double computeNumericalGradHy(const double x,
return gradHy;
}
+// Computes the 'gx' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradGx' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradGx'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
double computeNumericalGradGx(const double x,
const double y,
const elementStruct &element,
@@ -299,6 +310,9 @@ double computeNumericalGradGx(const double x,
return gradGx;
}
+// Computes the 'gy' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradGy' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradGy'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
double computeNumericalGradGy(const double x,
const double y,
const elementStruct &element,
diff --git a/srcs/BEM/functionsBEM.hpp b/srcs/BEM/functionsBEM.hpp
index c2cce53..d15d925 100644
--- a/srcs/BEM/functionsBEM.hpp
+++ b/srcs/BEM/functionsBEM.hpp
@@ -1,24 +1,20 @@
#include <gmsh.h>
#include <iostream>
-#include <map>
#include <sstream>
+#include <map>
#include <Eigen/Dense>
-#include <cmath>
-#include <iomanip>
-#include <stdio.h>
-#include <chrono>
#include <omp.h>
// Each element has all its information gathered as in this structure.
struct elementStruct{
- 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 midNodeTag;
- double length;
- double x1, y1, x2, y2;
- double midX, midY;
- double deltaX, deltaY;
+ std::size_t tag; // tag of the element.
+ double bcValue [2]; // values of [phi, phi_n].
+ size_t nodeTags [2]; // tags of the [first, second] nodes of the element.
+ size_t midNodeTag; // tag of the mid element node (or 0 if 2-node element).
+ double length; // length of the element.
+ double x1, y1, x2, y2; // coordinates of the first/second nodes of the element.
+ double midX, midY; // mid point coordinates.
+ double deltaX, deltaY; // x and y variations (sign matters).
};
/*
@@ -65,16 +61,17 @@ 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).
+// current BEM domain, with the appropriate information for each element (the bcValue will be computed later). Also returns
+// "domainPhysicalGroupDim" and "domainPhysicalGroupTag", the dimension and the tag (respectivelly) of the BEM domain of
+// interest. "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,
+ int &domainPhysicalGroupDim,
+ int &domainPhysicalGroupTag,
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
+// 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.
@@ -110,45 +107,56 @@ void computeElectrostaticPressure(std::map<int, double> &electrostaticPressure,
* Post processing
*/
-int getModel(std::vector<size_t> &domainNodeTags,
- std::map<std::size_t, std::size_t> &domainNodeIndices,
- std::vector<double> &domainNodeCoords,
- std::vector<int> &domainEntityTags,
- std::map<int, std::size_t> &domainEntityIndices,
- std::vector<std::vector<int>> &domainElementTypes,
- std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
- std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
- std::map<std::string, std::pair<int, int>> &physicalGroups,
- const std::string domainName);
-
-void computeInternalPhi(std::map<int, double> &internalNodalPhi,
- const std::vector<double> &domainNodeCoords,
- const std::vector<size_t> &domainNodeTags,
- std::map<std::size_t, std::size_t> &domainNodeIndices,
- const std::vector<elementStruct> &elementVector,
- const std::string integrationType,
- std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
- const bool &solutionType);
+// Retrieves every usefull information from the inside domain for the post-processing. "domainNodeTags" and "domainNodeCoords"
+// are vectors containing all the node tags and coordinates (respectivelly) of the domain. "domainNodeIndices" is a map
+// associating the index of each node tag for accessing the "domainNodeCoords" vector. "domainEntityTags" is a vector containing
+// the tag of each entity of the domain. "domainElementTypes" stores the types of the elements of the the domain. The first
+// dimension of this vector represents the different entities and the second, the different element types. "domainElementTags"
+// stores the tags of the elements of the domain. The first dimension of this vector represents the different entities, the
+// second, the different element types and the third, the different element tags. The "domainElementNodeTags" map associates to
+// each element tag a vector of node tags representing the nodes of that element. "domainPhysicalGroupTag" and
+// "domainPhysicalGroupDim" are the dimension and tag (respectivelly) of the BEM domain of interrest.
+void 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::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,
+ const int domainPhysicalGroupDim,
+ const int domainPhysicalGroupTag);
+// From the "elementVector" vector containing the information of all elements, constructs the "boundaryNodeToElements" map
+// associating to each boundary node tag, a vector of element tags of the contiguous elements. The vector of element tags
+// contains 2 element excepts in case of a mid element node (for higher order meshes), where only one element tag is stored.
void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
const std::vector<elementStruct> &elementVector);
-void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
- const std::vector<std::vector<int>> &domainElementTypes,
- const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
- const std::vector<elementStruct> &elementVector,
- const int domainPhysicalGroupTag);
-
-void phiInitialization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
- std::map<std::size_t, std::vector<std::size_t>> &position,
- const std::vector<std::vector<int>> &domainElementTypes,
- const std::vector<std::vector<std::vector<size_t>>> &domainElementTags);
-
-void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
- const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
- std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
- std::map<int, double> &internalNodalPhi);
+// Computes and stores the potential of each interior node in the "internalNodalPhi" (it can be accessed via the tag of the
+// node). "domainNodeTags" and "domainNodeCoords" are vectors containing all the node tags and coordinates (respectivelly) of the
+// domain. "domainNodeIndices" is a map associating the index of each node tag for accessing the "domainNodeCoords" vector.
+// "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM
+// domain. "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).
+// "boundaryNodeToElements" is a map associating to each boundary node tag, a vector of element tags of the contiguous elements.
+void computeInternalPhi(std::map<int, double> &internalNodalPhi,
+ const std::vector<double> &domainNodeCoords,
+ const std::vector<std::size_t> &domainNodeTags,
+ std::map<std::size_t, std::size_t> &domainNodeIndices,
+ const std::vector<elementStruct> &elementVector,
+ const std::string integrationType,
+ std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
+ const bool &solutionType);
+// Binds the potential at each node to the right elements. "phi" is a data structure storing the potiential of each node (fourth
+// dimension) of each element tag (third dimension) of each element type (second dimension) of each entity tag (first dimension).
+// "domainElementTags" stores the tag of each element (third dimension) of each element type (second dimension) of each entity
+// tag (first dimension). "domainElementTypes" stores the type of each element (second dimension) of each entity (first
+// dimension). The "domainElementNodeTags" map associates to each element tag a vector of node tags representing the nodes of
+// that element. "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of
+// the current BEM domain. "internalNodalPhi" is a map associating the value of the potential to each node tag of the domain
+// (boundaries excuded). "domainPhysicalGroupTag" is the tag of the physical domain of interrest. "boundaryNodeToElements" is a
+// map associating to each boundary node tag, a vector of element tags of the contiguous elements.
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,
@@ -158,6 +166,51 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
const int domainPhysicalGroupTag,
std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements);
+// Initializes the "phi" data structure in order for it to have the proper dimensions. It is a data structure storing the
+// potiential of each node (fourth dimension) of each element tag (third dimension) of each element type (second dimension) of
+// each entity tag (first dimension). The "position" map associates to each of the element tag of the domain the 3 indices for
+// accessing it in the "domainElementTags" or "phi" data structure. "domainElementTags" stores the tag of each element (third
+// dimension) of each element type (second dimension) of each entity tag (first dimension). "domainElementTypes" stores the type
+// of each element (second dimension) of each entity (first dimension).
+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);
+
+// Fills the "domainBoundaryElementTags" vector containing the tag of every elements possessing at least one node in common with
+// the boundary. "domainElementTags" stores the tag of each element (third dimension) of each element type (second dimension) of
+// each entity tag (first dimension). "domainElementTypes" stores the type of each element (second dimension) of each entity
+// (first dimension). "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file)
+// of the current BEM domain. "domainPhysicalGroupTag" is the tag of the physical domain of interrest.
+void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
+ const std::vector<std::vector<int>> &domainElementTypes,
+ const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
+ const std::vector<elementStruct> &elementVector,
+ const int domainPhysicalGroupTag);
+
+// Associates the potential to the nodes of each elements (if not already done) in the "phi" data structure. It is a data
+// structure storing the potiential of each node (fourth dimension) of each element tag (third dimension) of each element type
+// (second dimension) of each entity tag (first dimension). "domainElementTags" stores the tag of each element (third dimension)
+// of each element type (second dimension) of each entity tag (first dimension). The "domainElementNodeTags" map associates to
+// each element tag a vector of node tags representing the nodes of that element. "internalNodalPhi" is a map associating the
+// value of the potential to each node tag of the domain (boundaries excuded).
+void elementNodalPhiAssociation(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);
+
+
+// Computes the potential (if "computePhi" is set to 'true') and the electric field (if "computeElectricField" is set to 'true')
+// and stores them in the "phi" and "electricField" datastructures respectivelly. These assign the value of the potential /
+// electric field to each element tag (third dimension) of each element type (second dimension) of each entity tag (first
+// dimension) of the domain. The fourth dimension is either of size 1 (for the potential) or 3 (for the electric field).
+// "domainEntityTags" is a vector containing the tag of each entity of the domain. "domainElementTypes" stores the types of the
+// elements of the the domain. The first dimension of this vector represents the different entities and the second, the different
+// element types. "domainElementTags" stores the tags of the elements of the domain. The first dimension of this vector
+// represents the different entities, the second, the different element types and the third, the different element tags.
+// "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM
+// domain. "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 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,
@@ -169,11 +222,17 @@ void computeElementData(std::vector<std::vector<std::vector<std::vector<double>>
const bool computePhi,
const bool computeElectricField);
+// Deals with the display in gmsh. "viewTag" is the view tag of the data. "modelName" is the name of the model. "fileName" is the
+// name of the output file. "dataType" is either "ElementData" or "ElementNodeData" (which are the supported type of this
+// application). "domainElementTags" stores the tags of the elements of the domain. The first dimension of this vector represents
+// the different entities, the second, the different element types and the third, the different element tags. "data" is a data
+// structure storing the data to be displayed. The first three dimensions follow the same structure as the one of the
+// "elementNodeTags" and the fourth one is contains what needs to be displayed. "viewIndex" if the index of the view.
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::size_t>>> &domainElementTags,
const std::vector<std::vector<std::vector<std::vector<double>>>> &data,
const int viewIndex);
@@ -181,29 +240,101 @@ void displayData(const int viewTag,
* Verification
*/
+// Displays the time taken between the values of "start" and "end" (in seconds) if "showTime" is set to 'true'. "functionName" is
+// the name of the function that has been timed.
void displayTime(const double start,
- const double end,
- std::string functionName, const bool showTime);
+ const double end,
+ std::string functionName, const bool showTime);
+// Displays the different errors induced by the boundary element method compared with a simple analytical solution on a certain
+// model geometry ("BEM_convergence.geo"). "domainEntityTags" is a vector containing the tag of each entity of the domain.
+// "domainElementTypes" stores the types of the elements of the the domain. The first dimension of this vector represents the
+// different entities and the second, the different element types. "domainElementTags" stores the tags of the elements of the
+// domain. The first dimension of this vector represents the different entities, the second, the different element types and the
+// third, the different element tags. "elementalPhi" and "elementalElectricField" are the used data for comparing numerical and
+// analytical solutions. The convergence function can only be used with the "ElementData" data type.
void convergenceFun(const std::vector<int> &domainEntityTags,
const std::vector<std::vector<int>> &domainElementTypes,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalPhi,
const std::vector<std::vector<std::vector<std::vector<double>>>> &elementalElectricField);
-
/*
* Computations
*/
-double computeAnalyticalH(const double x, const double y, const elementStruct &element);
-double computeNumericalG(const double x, const double y, const elementStruct &element, const std::string integrationType);
-double computeNumericalGradHx(const double x, const double y, const elementStruct &element, const std::string integrationType);
-double computeNumericalGradHy(const double x, const double y, const elementStruct &element, const std::string integrationType);
-double computeNumericalGradGx(const double x, const double y, const elementStruct &element, const std::string integrationType);
-double computeNumericalGradGy(const double x, const double y, const elementStruct &element, const std::string integrationType);
-double computeAnalyticalGradHx(const double x, const double y, const elementStruct &element);
-double computeAnalyticalGradHy(const double x, const double y, const elementStruct &element);
-double computeAnalyticalGradGx(const double x, const double y, const elementStruct &element);
-double computeAnalyticalGradGy(const double x, const double y, const elementStruct &element);
-double computeAnalyticalG(const double x, const double y, const elementStruct &element);
+// Computes the 'h' element analytically. "x" and "y" are the coordinates at which 'h' should be computed. "element" is the
+// element 'j' while will helps to compute 'h'.
+double computeAnalyticalH(const double x,
+ const double y,
+ const elementStruct &element);
+
+// Computes the 'g' element analytically. "x" and "y" are the coordinates at which 'g' should be computed. "element" is the
+// element 'j' while will helps to compute 'g'.
+double computeAnalyticalG(const double x,
+ const double y,
+ const elementStruct &element);
+
+// Computes the 'hx' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradHx' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradHx'.
+double computeAnalyticalGradHx(const double x,
+ const double y,
+ const elementStruct &element);
+
+// Computes the 'hy' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradHy' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradHy'.
+double computeAnalyticalGradHy(const double x,
+ const double y,
+ const elementStruct &element);
+
+// Computes the 'gx' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradGx' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradGx'.
+double computeAnalyticalGradGx(const double x,
+ const double y,
+ const elementStruct &element);
+
+// Computes the 'gy' component of the gradient analytically. "x" and "y" are the coordinates at which 'gradGy' should be
+// computed. "element" is the element 'j' while will helps to compute 'gradGy'.
+double computeAnalyticalGradGy(const double x,
+ const double y,
+ const elementStruct &element);
+
+// Computes the 'g' element numerically. "x" and "y" are the coordinates at which 'g' should be computed. "element" is the
+// element 'j' while will helps to compute 'g'. "integrationType" is the type of numerical integration used (in case of numerical
+// computation).
+double computeNumericalG(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType);
+
+// Computes the 'hx' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradHx' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradHx'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
+double computeNumericalGradHx(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType);
+
+// Computes the 'hy' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradHy' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradHy'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
+double computeNumericalGradHy(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType);
+
+// Computes the 'gx' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradGx' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradGx'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
+double computeNumericalGradGx(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType);
+
+// Computes the 'gy' component of the gradient numerically. "x" and "y" are the coordinates at which 'gradGy' should be computed.
+// "element" is the element 'j' while will helps to compute 'gradGy'. "integrationType" is the type of numerical integration used
+// (in case of numerical computation).
+double computeNumericalGradGy(const double x,
+ const double y,
+ const elementStruct &element,
+ const std::string integrationType);
\ No newline at end of file
diff --git a/srcs/BEM/mainBEM.cpp b/srcs/BEM/mainBEM.cpp
index 0203578..cd639bb 100644
--- a/srcs/BEM/mainBEM.cpp
+++ b/srcs/BEM/mainBEM.cpp
@@ -11,19 +11,21 @@ int main(int argc, char **argv)
return 0;
}
+ int nthreads = omp_get_max_threads();
+ omp_set_num_threads(nthreads);
+
gmsh::initialize();
gmsh::open(argv[1]);
gmsh::model::mesh::generate(2);
- // 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(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, 1, meshUntangler);
if(postProcessing)
gmsh::fltk::run();
diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp
index c585455..691daac 100644
--- a/srcs/BEM/postProcessing.cpp
+++ b/srcs/BEM/postProcessing.cpp
@@ -1,46 +1,43 @@
-#include <gmsh.h>
-#include <iostream>
-#include <map>
-#include <sstream>
-#include <Eigen/Dense>
-#include <cmath>
-#include <iomanip>
-#include <stdio.h>
-#include <chrono>
#include "functionsBEM.hpp"
-#include <omp.h>
-#include <algorithm>
-#include <list>
-
-int getModel(std::vector<std::size_t> &domainNodeTags,
- std::map<std::size_t, std::size_t> &domainNodeIndices,
- std::vector<double> &domainNodeCoords,
- std::vector<int> &domainEntityTags,
- std::map<int, std::size_t> &domainEntityIndices,
- std::vector<std::vector<int>> &domainElementTypes,
- std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
- std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
- std::map<std::string, std::pair<int, int>> &physicalGroups,
- const std::string domainName)
+
+// Retrieves every usefull information from the inside domain for the post-processing. "domainNodeTags" and "domainNodeCoords"
+// are vectors containing all the node tags and coordinates (respectivelly) of the domain. "domainNodeIndices" is a map
+// associating the index of each node tag for accessing the "domainNodeCoords" vector. "domainEntityTags" is a vector containing
+// the tag of each entity of the domain. "domainElementTypes" stores the types of the elements of the the domain. The first
+// dimension of this vector represents the different entities and the second, the different element types. "domainElementTags"
+// stores the tags of the elements of the domain. The first dimension of this vector represents the different entities, the
+// second, the different element types and the third, the different element tags. The "domainElementNodeTags" map associates to
+// each element tag a vector of node tags representing the nodes of that element. "domainPhysicalGroupTag" and
+// "domainPhysicalGroupDim" are the dimension and tag (respectivelly) of the BEM domain of interrest.
+void 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::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,
+ const int domainPhysicalGroupDim,
+ const int domainPhysicalGroupTag)
{
- // Retrieves everything only for the domain of interrest.
- int physicalGroupDim = physicalGroups[domainName].first;
- int physicalGroupTag = physicalGroups[domainName].second;
+ // Only considers the domain of interrest.
+ int physicalGroupDim = domainPhysicalGroupDim;
+ int physicalGroupTag = domainPhysicalGroupTag;
+ // Retrieves the node tags and coordinates.
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.
+ // Constructs the "domainNodeIndices" map associating the index of each node tag in the "domainNodeCoords" vector.
for(std::size_t index = 0; index < domainNodeTags.size(); index++)
domainNodeIndices[domainNodeTags[index]] = index;
+ // Loops over the entities of the domain to fill the "elementTypes" and "elementTags" vectors and the "domainEntityIndices"
+ // and "domainElementNodeTags" maps.
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;
-
+ {
+ // Retrieves every 'element related' data from the tag of its entity.
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
std::vector<std::vector<std::size_t>> nodeTags;
@@ -50,8 +47,10 @@ int getModel(std::vector<std::size_t> &domainNodeTags,
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);
+ // Initializes a vector that will contain the tags of the nodes of the current element (the number of nodes only
+ // depends on the element type).
+ int nbNodesPerElement = nodeTags[k].size() / elementTags[k].size();
+ std::vector<std::size_t> elementNodeTags(nbNodesPerElement);
domainElementTypes[j][k] = elementTypes[k];
domainElementTags[j][k].resize(elementTags[k].size());
@@ -59,25 +58,29 @@ int getModel(std::vector<std::size_t> &domainNodeTags,
for(std::size_t l = 0; l < elementTags[k].size(); l++)
{
domainElementTags[j][k][l] = elementTags[k][l];
+
+ // Find the proper node tags of the elements.
int n = 0;
- for(std::size_t m = l*nbNodes; m < (l+1)*nbNodes; m++)
+ for(std::size_t m = l*nbNodesPerElement; m < (l+1)*nbNodesPerElement; m++)
{
elementNodeTags[n] = nodeTags[k][m];
n++;
}
+ // Stores it in the "domainElementNodeTags" map (element tag -> vector of node tags).
domainElementNodeTags[elementTags[k][l]] = elementNodeTags;
}
}
}
-
- return physicalGroupTag;
}
-
+// From the "elementVector" vector containing the information of all elements, constructs the "boundaryNodeToElements" map
+// associating to each boundary node tag, a vector of element tags of the contiguous elements. The vector of element tags
+// contains 2 element excepts in case of a mid element node (for higher order meshes), where only one element tag is stored.
void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements,
const std::vector<elementStruct> &elementVector)
{
+ // Loops over the elements and add the element tag to each of its corresponding nodes.
for(std::size_t i = 0; i < elementVector.size(); i++)
{
if(!boundaryNodeToElements.count(elementVector[i].nodeTags[0]))
@@ -95,6 +98,13 @@ void computeNodeToElements(std::map<std::size_t, std::vector<std::size_t>> &boun
}
}
+// Computes and stores the potential of each interior node in the "internalNodalPhi" (it can be accessed via the tag of the
+// node). "domainNodeTags" and "domainNodeCoords" are vectors containing all the node tags and coordinates (respectivelly) of the
+// domain. "domainNodeIndices" is a map associating the index of each node tag for accessing the "domainNodeCoords" vector.
+// "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM
+// domain. "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).
+// "boundaryNodeToElements" is a map associating to each boundary node tag, a vector of element tags of the contiguous elements.
void computeInternalPhi(std::map<int, double> &internalNodalPhi,
const std::vector<double> &domainNodeCoords,
const std::vector<std::size_t> &domainNodeTags,
@@ -109,18 +119,22 @@ void computeInternalPhi(std::map<int, double> &internalNodalPhi,
double phi = 0;
double h = 0;
double g = 0;
- std::vector<double> nodeVectorDomain(domainNodeTags.size(), 0);
+ std::vector<double> internalNodalPhiTmp(domainNodeTags.size(), 0); // Used for parallelization.
+ // Loops over all the internal nodes coordinates, computes the potential at that point and stores it in
+ // "internalNodalPhiTmp".
#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++)
{
+ // Do not consider the boundary nodes.
if(!boundaryNodeToElements.count(domainNodeTags[i]))
{
- // Coordinates of the node
+ // Gets the coordinates of the node.
x = domainNodeCoords[3*domainNodeIndices[domainNodeTags[i]]];
y = domainNodeCoords[3*domainNodeIndices[domainNodeTags[i]] + 1];
phi = 0;
+ // Loop over all the boundary element in order to compute the potential.
for(size_t j = 0; j < elementVector.size(); j++)
{
h = computeAnalyticalH(x, y, elementVector[j]);
@@ -131,14 +145,25 @@ void computeInternalPhi(std::map<int, double> &internalNodalPhi,
phi += h*elementVector[j].bcValue[0] - g*elementVector[j].bcValue[1];
}
- nodeVectorDomain[i] = phi;
+ internalNodalPhiTmp[i] = phi;
}
}
+ // Finally writes the potential in the "internalNodePhi" map.
for(std::size_t i = 0; i < domainNodeTags.size(); i++)
- internalNodalPhi[domainNodeTags[i]] = nodeVectorDomain[i];
+ internalNodalPhi[domainNodeTags[i]] = internalNodalPhiTmp[i];
}
+
+// Binds the potential at each node to the right elements. "phi" is a data structure storing the potiential of each node (fourth
+// dimension) of each element tag (third dimension) of each element type (second dimension) of each entity tag (first dimension).
+// "domainElementTags" stores the tag of each element (third dimension) of each element type (second dimension) of each entity
+// tag (first dimension). "domainElementTypes" stores the type of each element (second dimension) of each entity (first
+// dimension). The "domainElementNodeTags" map associates to each element tag a vector of node tags representing the nodes of
+// that element. "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of
+// the current BEM domain. "internalNodalPhi" is a map associating the value of the potential to each node tag of the domain
+// (boundaries excuded). "domainPhysicalGroupTag" is the tag of the physical domain of interrest. "boundaryNodeToElements" is a
+// map associating to each boundary node tag, a vector of element tags of the contiguous elements.
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,
@@ -147,22 +172,19 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
std::map<int, double> &internalNodalPhi,
const int domainPhysicalGroupTag,
std::map<std::size_t, std::vector<std::size_t>> &boundaryNodeToElements)
-{
- // This map may be replaced in the future. It allows to retain the position of the elements of the elementTag2d loop and linked their position to the loop of the elementTag2d_boundary (the tags of the elements in common with the boundary).
- // This map is the link between the array of tags of the 2d elements of the boundary and that of the tags of all the 2d elements.
+{
std::map<std::size_t, std::vector<std::size_t>> position;
phiInitialization(phi, position, domainElementTypes, domainElementTags);
- // // The first part of the algorithm consists in recovering the tag of the 2d elements of the domain which have at least one node in common with the boundary
- // // We initialise the first dimension of the vector that will contain all 2d elements in common with the boundary with the number of element types present in the domain.
std::vector<size_t> domainBoundaryElementTags;
boundaryElementTags(domainBoundaryElementTags, domainElementTypes, domainElementTags, elementVector, domainPhysicalGroupTag);
- // The second part of the algorithm will consist in assigning the correct potential value to the nodes of the elements in common with the boundary.
- // The first loop is on the elementType in 2d in common with the boundary.
+ // Loops over all the element that have node(s) in common with the boundary, and computes (or decides) what is the value of
+ // the potential at these boundary nodes.
#pragma omp parallel for schedule(guided)
for(std::size_t i = 0; i < domainBoundaryElementTags.size(); i++)
{
+ // Access everything related to the current element.
std::size_t elementTag = domainBoundaryElementTags[i];
int entityTagIndex = position[elementTag][0];
int elementTypeIndex = position[elementTag][1];
@@ -170,56 +192,50 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
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.
+ // Loops on the nodes of each boundary element.
for(int j = 0; j < nbNodes; j++)
{
- std::size_t elementNodeTag = nodeTags[j];
+ std::size_t elementNodeTag = nodeTags[j]; // Tag of each node.
+ // Only interrested in boundary nodes.
if(boundaryNodeToElements.count(elementNodeTag))
{
+ // Retrieves the tags of the contiguous boundary elements of the node.
std::vector<std::size_t> commonBoundaryElementIndices = boundaryNodeToElements[elementNodeTag];
- // Mid node of a segment on the boundary => take the value of the element.
+ // Mid node of a boundary element => take the value of the element.
if(commonBoundaryElementIndices.size() == 1)
- {
phi[entityTagIndex][elementTypeIndex][elementTagIndex][j] = elementVector[commonBoundaryElementIndices[0]].bcValue[0];
- }
- else if(commonBoundaryElementIndices.size() == 2)
+ // Several cases that lead to either select the value of the potential on one element or computes the average of
+ // the potential of two adjacent elements.
+ else
{
double currentPhi = 0.0;
int count = 0;
- // First element.
+ // Check the first contiguous 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 "adjacentNodeTag" is in "nodeTags"
if(std::count(nodeTags.begin(), nodeTags.end(), adjacentNodeTag))
- {
count++;
- }
+
currentPhi += elementVector[adjacentElementIndex].bcValue[0];
-
- // Second element.
+ // Check the second contiguous element.
adjacentElementIndex = commonBoundaryElementIndices[1];
adjacentNodeTag = elementVector[adjacentElementIndex].nodeTags[0];
if(adjacentNodeTag == elementNodeTag)
- {
adjacentNodeTag = elementVector[adjacentElementIndex].nodeTags[1];
- }
- // if 'adjacentNodeTag' is in 'nodeTags'
+ // if "adjacentNodeTag" is in "nodeTags"
if(std::count(nodeTags.begin(), nodeTags.end(), adjacentNodeTag))
{
if(count == 0)
- {
currentPhi = elementVector[adjacentElementIndex].bcValue[0];
- }
else
{
currentPhi += elementVector[adjacentElementIndex].bcValue[0];
@@ -227,13 +243,11 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
}
}
else
- {
if(count == 0)
{
currentPhi += elementVector[adjacentElementIndex].bcValue[0];
currentPhi /= 2;
}
- }
phi[entityTagIndex][elementTypeIndex][elementTagIndex][j] = currentPhi;
}
@@ -241,98 +255,126 @@ void elementNodeData(std::vector<std::vector<std::vector<std::vector<double>>>>
}
}
- //The last step is to loop over all the 2d elements in the domain. If the HUGE_VAL value is still present, this means that the element node does not belong to the boundary (since the value of the potential has not been modified) and is therefore part of the domain. We then use the map which associates the value of the potential with the tag of the node inside the domain.
- boundaryElementFinalization(phi, domainElementTags, domainElementNodeTags, internalNodalPhi);
+ elementNodalPhiAssociation(phi, domainElementTags, domainElementNodeTags, internalNodalPhi);
}
-
+// Initializes the "phi" data structure in order for it to have the proper dimensions. It is a data structure storing the
+// potiential of each node (fourth dimension) of each element tag (third dimension) of each element type (second dimension) of
+// each entity tag (first dimension). The "position" map associates to each of the element tag of the domain the 3 indices for
+// accessing it in the "domainElementTags" or "phi" data structure. "domainElementTags" stores the tag of each element (third
+// dimension) of each element type (second dimension) of each entity tag (first dimension). "domainElementTypes" stores the type
+// of each element (second dimension) of each entity (first dimension).
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)
{
+ // First dimension: entity tags.
phi.resize(domainElementTypes.size());
for(std::size_t i = 0; i < domainElementTypes.size(); i++)
{
+ // Second dimension: element types.
phi[i].resize(domainElementTypes[i].size());
#pragma omp parallel for schedule(guided) shared(phi, domainElementTypes, domainElementTags, position)
for(std::size_t j = 0; j < domainElementTypes[i].size(); j++)
{
+ // "getElementProperties" is called in order to access the number of nodes of each element (only depends on the
+ // element type).
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);
+ gmsh::model::mesh::getElementProperties(domainElementTypes[i][j], elementName, dim, order, nbNodes, localNodeCoord,
+ nbPrimaryNodes);
+ // Third dimension: element tags.
phi[i][j].resize(domainElementTags[i][j].size());
for(std::size_t k = 0; k < domainElementTags[i][j].size(); k++)
{
- phi[i][j][k].resize(nbNodes, HUGE_VAL);
+ // Fourth dimension: node tags.
+ phi[i][j][k].resize(nbNodes, HUGE_VAL); // First set to a high value
std::vector<std::size_t> pos = {i, j, k};
#pragma omp critical
- position[domainElementTags[i][j][k]] = pos;
+ position[domainElementTags[i][j][k]] = pos; // Stores the indices
}
}
}
}
+// Fills the "domainBoundaryElementTags" vector containing the tag of every elements possessing at least one node in common with
+// the boundary. "domainElementTags" stores the tag of each element (third dimension) of each element type (second dimension) of
+// each entity tag (first dimension). "domainElementTypes" stores the type of each element (second dimension) of each entity
+// (first dimension). "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file)
+// of the current BEM domain. "domainPhysicalGroupTag" is the tag of the physical domain of interrest.
void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
const std::vector<std::vector<int>> &domainElementTypes,
const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
const std::vector<elementStruct> &elementVector,
const int domainPhysicalGroupTag)
{
- #pragma omp parallel
+ // #pragma omp parallel
{
- std::vector<std::size_t> localElementTags;
+ std::vector<std::size_t> localElementTags; // For parallelization.
- // 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)
+ // Loops over each boundary element and retrieves the element in contact with its boundary nodes.
+ // #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.
+ // Adds the elements tags to the local vector.
for(std::size_t j = 0; j < elementTags.size(); j++)
{
+ // Retrieves the physical group to which the element belongs because it is possible that the element belongs to
+ // an adjacent physical surface.
int elementType;
std::vector<size_t> nodeTags;
int entityDim;
int entityTag;
std::vector<int> physicalTags;
- #pragma omp critical
+ // #pragma omp critical
{
gmsh::model::mesh::getElement(elementTags[j], elementType, nodeTags, entityDim, entityTag);
gmsh::model::getPhysicalGroupsForEntity(entityDim, entityTag, physicalTags);
}
+ // Add the element tag (to a local 'thread safe' vector) if the physical tag is the right one.
for(std::size_t k = 0; k < physicalTags.size(); k++)
if(physicalTags[k] == domainPhysicalGroupTag)
localElementTags.push_back(elementTags[j]);
}
}
- #pragma omp critical
+ // // Fills the "domainBoundaryElementTags" vector.
+ // #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.
+ // Allows duplicate tags to be removed. If an element has several nodes in common with the boundary, it has been counted
+ // several times before.
sort( domainBoundaryElementTags.begin(), domainBoundaryElementTags.end() );
- domainBoundaryElementTags.erase( unique( domainBoundaryElementTags.begin(), domainBoundaryElementTags.end() ), domainBoundaryElementTags.end() );
+ domainBoundaryElementTags.erase( unique( domainBoundaryElementTags.begin(), domainBoundaryElementTags.end() ),
+ domainBoundaryElementTags.end() );
}
-
-void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector<double>>>> &phi,
- const std::vector<std::vector<std::vector<std::size_t>>> &domainElementTags,
- std::map<std::size_t, std::vector<std::size_t>> &domainElementNodeTags,
- std::map<int, double> &internalNodalPhi)
+// Associates the potential to the nodes of each elements (if not already done) in the "phi" data structure. It is a data
+// structure storing the potiential of each node (fourth dimension) of each element tag (third dimension) of each element type
+// (second dimension) of each entity tag (first dimension). "domainElementTags" stores the tag of each element (third dimension)
+// of each element type (second dimension) of each entity tag (first dimension). The "domainElementNodeTags" map associates to
+// each element tag a vector of node tags representing the nodes of that element. "internalNodalPhi" is a map associating the
+// value of the potential to each node tag of the domain (boundaries excuded).
+void elementNodalPhiAssociation(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)
{
+ // Loops over each nodes of the domain (entity tags, then element types, than element tags, the node tags).
for(std::size_t i = 0; i < domainElementTags.size(); i++)
{
for(std::size_t j = 0; j < domainElementTags[i].size(); j++)
@@ -344,6 +386,9 @@ 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 the potential has not already been modified in "phi" and if the node does not belong to the boundary,
+ // retrieves the potential of that node thanks to the "interNodalPhi".
if(phi[i][j][k][l] == HUGE_VAL && internalNodalPhi.count(elementNodeTag))
phi[i][j][k][l] = internalNodalPhi[elementNodeTag];
}
@@ -352,7 +397,17 @@ void boundaryElementFinalization(std::vector<std::vector<std::vector<std::vector
}
}
-
+// Computes the potential (if "computePhi" is set to 'true') and the electric field (if "computeElectricField" is set to 'true')
+// and stores them in the "phi" and "electricField" datastructures respectivelly. These assign the value of the potential /
+// electric field to each element tag (third dimension) of each element type (second dimension) of each entity tag (first
+// dimension) of the domain. The fourth dimension is either of size 1 (for the potential) or 3 (for the electric field).
+// "domainEntityTags" is a vector containing the tag of each entity of the domain. "domainElementTypes" stores the types of the
+// elements of the the domain. The first dimension of this vector represents the different entities and the second, the different
+// element types. "domainElementTags" stores the tags of the elements of the domain. The first dimension of this vector
+// represents the different entities, the second, the different element types and the third, the different element tags.
+// "elementVector" is a vector storing the boundary elements (see "elementStruct" in the 'functions.hpp' file) of the current BEM
+// domain. "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 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,
@@ -366,17 +421,23 @@ void computeElementData(std::vector<std::vector<std::vector<std::vector<double>>
{
std::vector<double> barycenters;
+ // Loops over all the entity tags and element types.
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());
+
+ // Retrieves the barycenters of each elements.
gmsh::model::mesh::getBarycenters(domainElementTypes[i][j], domainEntityTags[i], false, false, barycenters);
+
+ // Loops over each barycenter and compute the potential and electric field at that point.
+ #pragma omp parallel for schedule(static)
for(size_t k = 0; k < domainElementTags[i][j].size(); k++)
{
double x = barycenters[3*k];
@@ -386,58 +447,54 @@ void computeElementData(std::vector<std::vector<std::vector<std::vector<double>>
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;
+
+ // Loops over all the boundary elements in order to compute the potential and the electric field.
+ for(std::size_t el = 0; el < elementVector.size(); el++)
{
- double localCurrentPhi = 0;
- double localGradPhiX = 0;
- double localGradPhiY = 0;
+ if(computePhi)
+ {
+ h = computeAnalyticalH(x, y, elementVector[el]);
+ if(!solutionType)
+ g = computeNumericalG(x, y, elementVector[el], integrationType);
+ else
+ g = computeAnalyticalG(x, y, elementVector[el]);
- double h, hX, hY;
- double g, gX, gY;
-
- #pragma omp for schedule(guided)
- for(std::size_t el = 0; el < elementVector.size(); el++)
+ localCurrentPhi += h*elementVector[el].bcValue[0] - g*elementVector[el].bcValue[1];
+ }
+ if(computeElectricField)
{
- if(computePhi)
+ if(!solutionType)
{
- 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];
+ 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);
}
- if(computeElectricField)
+ else
{
- 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];
+ hX = computeAnalyticalGradHx(x, y, elementVector[el]);
+ gX = computeAnalyticalGradGx(x, y, elementVector[el]);
+ hY = computeAnalyticalGradHy(x, y, elementVector[el]);
+ gY = computeAnalyticalGradGy(x, y, elementVector[el]);
}
- }
- #pragma omp critical
- {
- currentPhi += localCurrentPhi;
- gradPhiX += localGradPhiX;
- gradPhiY += localGradPhiY;
+ localGradPhiX += hX*elementVector[el].bcValue[0] - gX*elementVector[el].bcValue[1];
+ localGradPhiY += hY*elementVector[el].bcValue[0] - gY*elementVector[el].bcValue[1];
}
}
+ // Updates the potential and electric field (for the parallelization)
+ currentPhi += localCurrentPhi;
+ gradPhiX += localGradPhiX;
+ gradPhiY += localGradPhiY;
+
+ // Writes the data in the data structures.
phi[i][j][k].resize(1, currentPhi);
electricField[i][j][k].resize(3);
electricField[i][j][k][0] = -gradPhiX;
@@ -448,6 +505,12 @@ void computeElementData(std::vector<std::vector<std::vector<std::vector<double>>
}
}
+// Deals with the display in gmsh. "viewTag" is the view tag of the data. "modelName" is the name of the model. "fileName" is the
+// name of the output file. "dataType" is either "ElementData" or "ElementNodeData" (which are the supported type of this
+// application). "domainElementTags" stores the tags of the elements of the domain. The first dimension of this vector represents
+// the different entities, the second, the different element types and the third, the different element tags. "data" is a data
+// structure storing the data to be displayed. The first three dimensions follow the same structure as the one of the
+// "elementNodeTags" and the fourth one is contains what needs to be displayed. "viewIndex" if the index of the view.
void displayData(const int viewTag,
const std::string modelName,
const std::string fileName,
@@ -460,9 +523,10 @@ void displayData(const int viewTag,
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::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);
+ gmsh::option::setNumber("View[" + std::to_string(viewIndex) + "].Visible", 0); // preventing the superimposed plotting
}
\ No newline at end of file
diff --git a/srcs/BEM/solverBEM.cpp b/srcs/BEM/solverBEM.cpp
index 67ab21d..ef14457 100644
--- a/srcs/BEM/solverBEM.cpp
+++ b/srcs/BEM/solverBEM.cpp
@@ -1,15 +1,4 @@
#include "functionsBEM.hpp"
-#include <gmsh.h>
-#include <iostream>
-#include <map>
-#include <sstream>
-#include <Eigen/Dense>
-#include <cmath>
-#include <iomanip>
-#include <stdio.h>
-#include <chrono>
-#include <omp.h>
-#include <list>
// 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
@@ -100,7 +89,9 @@ void solverBEM(std::map<int, double> &electrostaticPressure,
// 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);
+ singleDomainSolverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, iteration,
+ "BEM_domain_" + domainNames[i], phiViewTag, electricFieldViewTag, computePhi,
+ computeElectricField);
if(iteration == 1)
std::cout << "--------------------[BEM SOLVER]--------------------\n\n";
@@ -135,8 +126,10 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
// Parameters:
bool showTime = false; // true: displays the different time measurements.
- bool convergence = false; // true: computes the error between the analytical solution and the numerical one on a simple configuration (/!\ use this with the BEM_convergence.geo model).
- bool visualisation = true; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of the potential (Remark: electric field is always displayed with 'ElementData' type).
+ bool convergence = false; // true: computes the error between the analytical solution and the numerical one on a simple
+ // configuration (/!\ use this with the BEM_convergence.geo model).
+ bool visualisation = true; // true: 'ElementNodeData' visualization of the potential, false: 'ElementData' visualization of
+ // the potential (Remark: electric field is always displayed with 'ElementData' type).
bool solutionType = true; // true: analytical computations, false: numerical computations.
// Used for displaying the time of different sections.
@@ -152,8 +145,10 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
// Fills the vector of boundary elements.
start = omp_get_wtime();
- std::map<std::string, std::pair<int, int>> physicalGroups;
- std::size_t nbDirichletElements = fillElementVector(elementVector, physicalGroups, nodalDisplacements, domainName);
+ int domainPhysicalGroupDim;
+ int domainPhysicalGroupTag;
+ std::size_t nbDirichletElements = fillElementVector(elementVector, domainPhysicalGroupDim, domainPhysicalGroupTag,
+ nodalDisplacements, domainName);
end = omp_get_wtime();
displayTime(start, end, "'fillElementVector'", showTime);
@@ -201,12 +196,12 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
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;
- int domainPhysicalGroupTag = getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainEntityIndices, domainElementTypes, domainElementTags, domainElementNodeTags, physicalGroups, domainName);
+ getModel(domainNodeTags, domainNodeIndices, domainNodeCoords, domainEntityTags, domainElementTypes, domainElementTags,
+ domainElementNodeTags, domainPhysicalGroupDim, domainPhysicalGroupTag);
end = omp_get_wtime();
displayTime(start, end, "'getModel'", showTime);
@@ -224,46 +219,54 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
// Visualisation by ElementNodeData
if(visualisation)
{
- // Computes a map that stores 'nodeTag' (of the boundary) -> vector of adjacent 'elementIndices' (because only 1 element for higher order elements).
+ // Computes a map that stores 'nodeTag' (of the boundary) -> vector of adjacent 'elementIndices' (because only 1
+ // element for higher order elements).
std::map<std::size_t, std::vector<std::size_t>> boundaryNodeToElements;
computeNodeToElements(boundaryNodeToElements, elementVector);
- // Computes the value of the potential at each node strictly included in the domain and stores it in the "internalNodalPhi" map
+ // Computes the value of the potential at each node strictly included in the domain and stores it in the
+ // "internalNodalPhi" map.
start = omp_get_wtime();
std::map<int, double> internalNodalPhi;
- computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector, integrationType, boundaryNodeToElements, solutionType);
+ computeInternalPhi(internalNodalPhi, domainNodeCoords, domainNodeTags, domainNodeIndices, elementVector,
+ integrationType, boundaryNodeToElements, solutionType);
end = omp_get_wtime();
displayTime(start, end, "'computeInternalPhi'", showTime);
// Associates the value of the potential to each node of an element.
start = omp_get_wtime();
- elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector, internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements);
+ elementNodeData(elementNodalPhi, domainElementTags, domainElementTypes, domainElementNodeTags, elementVector,
+ internalNodalPhi, domainPhysicalGroupTag, boundaryNodeToElements);
end = omp_get_wtime();
displayTime(start, end, "'elementNodeData'", showTime);
// Computes the electric field (with "ElementData" type).
start = omp_get_wtime();
- computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, false, computeElectricField);
+ computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags,
+ elementVector, solutionType, integrationType, false, computeElectricField);
end = omp_get_wtime();
displayTime(start, end, "'computeElementData' (only the electric field)", showTime);
// Displays what is needed.
if(computePhi)
{
- displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi, viewIndex);
+ displayData(phiViewTag, names[0], "results.msh", "ElementNodeData", domainElementTags, elementNodalPhi,
+ viewIndex);
viewIndex++;
}
if(computeElectricField)
- displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex);
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags,
+ elementalElectricField, viewIndex);
}
// Visualisation by ElementData
else
{
// Computes the potential and the electric field.
start = omp_get_wtime();
- computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags, elementVector, solutionType, integrationType, computePhi, computeElectricField);
+ computeElementData(elementalPhi, elementalElectricField, domainElementTags, domainElementTypes, domainEntityTags,
+ elementVector, solutionType, integrationType, computePhi, computeElectricField);
end = omp_get_wtime();
displayTime(start, end, "'computeElementData'", showTime);
@@ -274,7 +277,8 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
viewIndex++;
}
if(computeElectricField)
- displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags, elementalElectricField, viewIndex);
+ displayData(electricFieldViewTag, names[0], "results.msh", "ElementData", domainElementTags,
+ elementalElectricField, viewIndex);
// Study of the error on a precise model.
if(convergence)
@@ -282,9 +286,13 @@ void singleDomainSolverBEM(std::map<int, double> &electrostaticPressure,
std::cout << "\tCONVERGENCE:\n";
std::string convergenceModelName = "BEM_convergence";
if(names[0].compare(convergenceModelName) == 0)
- convergenceFun(domainEntityTags, domainElementTypes, domainElementTags, elementalPhi, elementalElectricField);
+ convergenceFun(domainEntityTags, domainElementTypes, domainElementTags, elementalPhi,
+ elementalElectricField);
else
- std::cout << "\tWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '" << convergenceModelName + ".geo" << "' model if you want to check for convergence.\n\n";
+ {
+ std::cout << "\tWarning: Cannot check convergence on model '" << names[0] << ".geo'. Use the '";
+ std::cout << convergenceModelName + ".geo" << "' model if you want to check for convergence.\n\n";
+ }
}
}
}
diff --git a/srcs/BEM/verification.cpp b/srcs/BEM/verification.cpp
index 48d7bbd..c3136c8 100644
--- a/srcs/BEM/verification.cpp
+++ b/srcs/BEM/verification.cpp
@@ -1,50 +1,26 @@
-#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>
+// Displays the time taken between the values of "start" and "end" (in seconds) if "showTime" is set to 'true'. "functionName" is
+// the name of the function that has been timed.
void displayTime(const double start,
- const double end,
- std::string functionName, const bool showTime)
+ const double end,
+ std::string functionName, const bool showTime)
{
- //displays the time of the desired code section
if(showTime)
{
double duration = end - start;
std::cout << functionName << "\n";
std::cout << "Elapsed time: " << 1000000*duration << " [µs]\tor " << duration << " [s]\n\n";
-
- // std::cout << std::endl;
- // std::cout << function_name << std::endl;
- // std::cout << "Elapsed time in nanoseconds: "
- // << std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count()
- // << " ns" << std::endl;
-
- // std::cout << "Elapsed time in microseconds: "
- // << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
- // << " µs" << std::endl;
-
- // std::cout << "Elapsed time in milliseconds: "
- // << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
- // << " ms" << std::endl;
-
- // std::cout << "Elapsed time in seconds: "
- // << std::chrono::duration_cast<std::chrono::seconds>(end - start).count()
- // << " sec";
- // std::cout << std::endl;
- // std::cout << std::endl;
}
}
+// Displays the different errors induced by the boundary element method compared with a simple analytical solution on a certain
+// model geometry ("BEM_convergence.geo"). "domainEntityTags" is a vector containing the tag of each entity of the domain.
+// "domainElementTypes" stores the types of the elements of the the domain. The first dimension of this vector represents the
+// different entities and the second, the different element types. "domainElementTags" stores the tags of the elements of the
+// domain. The first dimension of this vector represents the different entities, the second, the different element types and the
+// third, the different element tags. "elementalPhi" and "elementalElectricField" are the used data for comparing numerical and
+// analytical solutions. The convergence function can only be used with the "ElementData" data type.
void convergenceFun(const std::vector<int> &domainEntityTags,
const std::vector<std::vector<int>> &domainElementTypes,
const std::vector<std::vector<std::vector<size_t>>> &domainElementTags,
--
GitLab