From b22ee5a3a8ab01d062194a3cfb573166e8c6c1b1 Mon Sep 17 00:00:00 2001
From: Louis Denis <louis.denis@student.uliege.be>
Date: Fri, 13 May 2022 11:56:31 +0200
Subject: [PATCH] cleaning FEM code + comments
---
srcs/FEM/complex_validation.geo | 14 +-
srcs/FEM/functionsFEM.cpp | 1786 +++++++++++++-----------
srcs/FEM/functionsFEM.hpp | 283 +++-
srcs/FEM/large_rotation_validation.geo | 2 +-
srcs/FEM/mainFEM.cpp | 13 +-
srcs/FEM/solverFEM.cpp | 474 ++-----
6 files changed, 1399 insertions(+), 1173 deletions(-)
diff --git a/srcs/FEM/complex_validation.geo b/srcs/FEM/complex_validation.geo
index 2c7851a..c6497d7 100644
--- a/srcs/FEM/complex_validation.geo
+++ b/srcs/FEM/complex_validation.geo
@@ -5,12 +5,12 @@ Ly = 2;
nx = 250;
ny = 100;
-Point(1) = {0, 0, 0, 0.1}; //4th number is mesh density
-Point(2) = {Lx, 0, 0, 0.2};
-Point(3) = {Lx, Ly, 0, 0.2};
-Point(4) = {0, Ly, 0, 0.1};
-Point(5) = {Lx/2, 0, 0, 0.15};
-Point(6) = {Lx/2, Ly, 0, 0.15};
+Point(1) = {0, 0, 0, 0.05}; //4th number is mesh density
+Point(2) = {Lx, 0, 0, 0.1};
+Point(3) = {Lx, Ly, 0, 0.1};
+Point(4) = {0, Ly, 0, 0.05};
+Point(5) = {Lx/2, 0, 0, 0.075};
+Point(6) = {Lx/2, Ly, 0, 0.075};
Line(1) = {1, 5};
Line(2) = {2, 5};
Line(3) = {3, 2};
@@ -56,4 +56,4 @@ SetNumber("Boundary Conditions/fixed_node/uy",2);
Physical Curve("BEM_FEM_boundary", 10) = {1,2,3,4,5,6};
-SetNumber("Non_linear_solver", 1);
\ No newline at end of file
+SetNumber("Non_linear_solver", 0);
\ No newline at end of file
diff --git a/srcs/FEM/functionsFEM.cpp b/srcs/FEM/functionsFEM.cpp
index e6cb0bd..4f69cc0 100644
--- a/srcs/FEM/functionsFEM.cpp
+++ b/srcs/FEM/functionsFEM.cpp
@@ -1,25 +1,27 @@
#ifdef _MSC_VER
-#include <gmsh.h_cwrap> // gmsh main header
+#include <gmsh.h_cwrap>
#else
-#include <gmsh.h> // gmsh main header
+#include <gmsh.h>
#endif
#include "functionsFEM.hpp"
-#include <iostream> // for std::cout
+#include <iostream>
#include <map>
-#include <sstream> // for std::stringstream
-#include <Eigen/Dense> // Eigen library
-#include <Eigen/Sparse> // Eigen library for sparse matrices
+#include <sstream>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
#include <Eigen/SparseCholesky> // solving sparse linear systems
#include <Eigen/Core>
#include <cmath>
-#include <algorithm> // to sort and merge the FEM node vectors
+#include <algorithm>
#ifdef _OPENMP
#include <omp.h>
#endif
#define PI M_PI
+// returns the boolean non_linear_solver parameter based on what is set in the .geo file,
+// add SetNumber("Non_linear_solver", 0); in .geo file to use linear solver.
bool getNonLinearParameter()
{
std::vector<std::string> keys;
@@ -44,42 +46,190 @@ bool getNonLinearParameter()
}
}
-//Compute the matrix H
-Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu)
+// fills the "FEM_nodeTags" vector with the nodes in the FEM domain, fills the "FEM_nodeCoordsMap" with the coordinates (x,y)
+// corresponding to each nodeTag, looking at entity tags "tags" of dimension "dim" corresponding to the FEM domain.
+void RetrieveFEMNodeTagsAndCoords(std::vector<std::size_t> & FEM_nodeTags,
+ std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap,
+ int & dim, std::vector<int> & tags)
{
- Eigen::Matrix<double, 3, 3> H_matrix;
- H_matrix(0,0) = 1.;
- H_matrix(1,0) = nu;
- H_matrix(2,0) = 0.;
- H_matrix(0,1) = nu;
- H_matrix(1,1) = 1.;
- H_matrix(2,1) = 0.;
- H_matrix(0,2) = 0.;
- H_matrix(1,2) = 0.;
- H_matrix(2,2) = (1-nu)/2;
- H_matrix = E/(1-nu*nu)*H_matrix;
- return H_matrix;
+ //loop over the entities of the FEM domain to retrieve all nodeTags by taking care not to take some nodes twice
+ //using sort and set_union methods + retrieving the node coordinates and storing them in a map
+ std::vector<double> FEM_nodecoord;
+ std::vector<double> FEM_nodeparametricCoord;
+ std::vector<std::vector<std::size_t>> tmp_nodeTags(tags.size()); // will store the nodeTags associated to one given entity of the domain
+ gmsh::model::mesh::getNodes(tmp_nodeTags[0], FEM_nodecoord, FEM_nodeparametricCoord, dim, tags[0], true);
+ for(std::size_t j = 0; j < tmp_nodeTags[0].size(); j++)
+ {
+ FEM_nodeCoordsMap[tmp_nodeTags[0][j]] = std::pair<double, double>(FEM_nodecoord[3*j], FEM_nodecoord[3*j+1]);
+ }
+ std::sort(tmp_nodeTags[0].begin(), tmp_nodeTags[0].end()); // need to sort the vectors in order to merge them afterwards
+ for (std::size_t i = 1; i < tags.size(); i++)
+ { // in the case the domain contains multiple entities
+ gmsh::model::mesh::getNodes(tmp_nodeTags[i], FEM_nodecoord, FEM_nodeparametricCoord, dim, tags[i], true);
+ for(std::size_t j = 0; j < tmp_nodeTags[i].size(); j++)
+ {
+ FEM_nodeCoordsMap[tmp_nodeTags[i][j]] = std::pair<double, double>(FEM_nodecoord[3*j], FEM_nodecoord[3*j+1]);
+ }
+ std::sort(tmp_nodeTags[i].begin(), tmp_nodeTags[i].end());
+ std::set_union(tmp_nodeTags[i-1].begin(), tmp_nodeTags[i-1].end(),
+ tmp_nodeTags[i].begin(), tmp_nodeTags[i].end(), std::back_inserter(FEM_nodeTags)); // works only for 2 domains -> see little trick just below
+
+ }
+ if(tags.size() == 1){ // if the domain contains one single entity
+ FEM_nodeTags = tmp_nodeTags[0];
+ }
+ else{ // little trick for more than 2 domains
+ std::sort(FEM_nodeTags.begin(), FEM_nodeTags.end());
+ auto last = std::unique(FEM_nodeTags.begin(), FEM_nodeTags.end());
+ FEM_nodeTags.erase(last, FEM_nodeTags.end());
+ }
}
-void computeBmatrix(const int &j, const int &numNodes, const std::vector<double> & basisFunctionsGradient,
- const std::vector<double> & determinants,
- const Eigen::Matrix<double, 2, 2> & jacobinvtrans,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & B_matrix)
+// fills the "FEM_elementTags" vector with all element tags of the elements of dimension "dim" (in practice: 2)
+// in the entities associated to the entity tags "tags" (corresponding to the FEM domain).
+// associates an initialized elementData structure to each element tag through the "FEM_kinematicsMap" map,
+// using the "FEM_nodeCoordsMap" map gathering the coordinates of each node in the FEM domain.
+int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags)
{
- for (int l = 0; l < numNodes; ++l){ // looping over the number of shape functions
- B_matrix(0,2*l) = jacobinvtrans(0,0)*basisFunctionsGradient[3*l+numNodes*3*j]+jacobinvtrans(0,1)*basisFunctionsGradient[3*l+numNodes*3*j+1];
- B_matrix(1,2*l) = 0.;
- B_matrix(2,2*l) = jacobinvtrans(1,0)*basisFunctionsGradient[3*l+numNodes*3*j]+jacobinvtrans(1,1)*basisFunctionsGradient[3*l+numNodes*3*j+1];
- B_matrix(0,2*l+1) = 0.;
- B_matrix(1,2*l+1) = B_matrix(2,2*l);
- B_matrix(2,2*l+1) = B_matrix(0,2*l);
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> elnodeTags;
+
+ std::string elementName;
+ int elDim, order, numNodes, numPrimaryNodes;
+ std::vector<double> localNodeCoord;
+
+ int nbelems=0; // will store total number of elements in FEM domain
+ for(std::size_t i=0; i< tags.size(); ++i)
+ {
+ gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[i]);
+ for(std::size_t j=0; j< elementTags.size(); ++j)
+ {
+ nbelems+=elementTags[j].size();
+ // filling the FEM_elementTags vector
+ FEM_elementTags.insert(FEM_elementTags.end(), elementTags[j].begin(), elementTags[j].end());
+
+ // getting element properties, in particular: number of nodes
+ gmsh::model::mesh::getElementProperties(elementTypes[j],
+ elementName, elDim, order,
+ numNodes, localNodeCoord,
+ numPrimaryNodes);
+ for(std::size_t k = 0; k < elementTags[j].size(); k++)
+ {
+ elementData element;
+ std::vector<size_t> tmp_nodes(numNodes);
+ std::vector<double> tmp_x(numNodes);
+ std::vector<double> tmp_y(numNodes);
+ std::vector<double> tmp_u(numNodes);
+ std::vector<double> tmp_v(numNodes);
+ std::vector<double> tmp_ui(numNodes);
+ std::vector<double> tmp_vi(numNodes);
+ double tmp_xc = 0;
+ double tmp_yc = 0;
+ for(int l = 0; l < numNodes; l++)
+ {
+ tmp_nodes[l] = elnodeTags[j][numNodes*k + l];
+ tmp_x[l] = FEM_nodeCoordsMap[tmp_nodes[l]].first;
+ tmp_y[l] = FEM_nodeCoordsMap[tmp_nodes[l]].second;
+ tmp_u[l] = 0; // initializing the displacement at zero
+ tmp_v[l] = 0;
+ tmp_ui[l] = 0;
+ tmp_vi[l] = 0;
+ tmp_xc += tmp_x[l];
+ tmp_yc += tmp_y[l];
+ }
+ tmp_xc = tmp_xc / numNodes;
+ tmp_yc = tmp_yc / numNodes;
+
+ std::vector<double> tmp_xi(numNodes);
+ std::vector<double> tmp_yi(numNodes);
+ for(int l = 0; l < numNodes; l++)
+ {
+ tmp_xi[l] = tmp_x[l] - tmp_xc;
+ tmp_yi[l] = tmp_y[l] - tmp_yc;
+ }
+ element.nodes = tmp_nodes;
+ element.x = tmp_x;
+ element.y = tmp_y;
+ element.xc = tmp_xc;
+ element.yc = tmp_yc;
+ element.xi = tmp_xi;
+ element.yi = tmp_yi;
+ element.u = tmp_u;
+ element.v = tmp_v;
+ element.uc = 0;
+ element.vc = 0;
+ element.theta = 0;
+ element.ui = tmp_ui;
+ element.vi = tmp_vi;
+ element.pl = Eigen::Matrix<double, Eigen::Dynamic, 1>::Zero(2*numNodes,1);
+ element.A = Eigen::Matrix<double, Eigen::Dynamic, 1>(2*numNodes,1);
+ element.G = Eigen::Matrix<double, 1, Eigen::Dynamic>(1,2*numNodes);
+ element.E = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(2*numNodes, 2*numNodes);
+
+ FEM_kinematicsMap[elementTags[j][k]] = element;
+ updateKinMatrices(& FEM_kinematicsMap[elementTags[j][k]], false);
+ }
+ }
+ }
+ return nbelems;
+}
+
+// updates the kinematics matrices of one element "element" based on its updated nodal positions and rotation angle.
+// if "computeFl" is set as true, also computes the local nodal forces vector "fl = Kl*pl", only OK if
+// the local stiffness matrix Kl has already been initialized, using "fillElementalKandF".
+void updateKinMatrices(elementData* element, bool computeFl)
+{
+ // for the different formulas, please refer to the report
+ size_t numNodes = element->nodes.size();
+ // update A matrix
+ for(size_t i = 0; i < numNodes; i++)
+ {
+ element->A(2*i, 0) = -(element->vi[i]+element->yi[i]);
+ element->A(2*i + 1, 0) = element->ui[i]+element->xi[i];
+ }
+
+ // update G matrix
+ double denominator = 0;
+ for(size_t i = 0; i < numNodes; i++)
+ {
+ element->G(0, 2*i) = -element->yi[i];
+ element->G(0, 2*i + 1) = element->xi[i];
+ denominator += element->xi[i]*(element->ui[i]+element->xi[i]) + element->yi[i]*(element->vi[i]+element->yi[i]);
+ }
+ element->G *= 1 / denominator;
+
+ // update P matrix
+ element->P = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(2*numNodes,2*numNodes) - element->A*element->G;
+
+ // update E matrix
+ double tmp_cos = cos(element->theta);
+ double tmp_sin = sin(element->theta);
+ for(size_t i = 0; i < numNodes; i++)
+ {
+ element->E(2*i, 2*i) = tmp_cos;
+ element->E(2*i, 2*i + 1) = -tmp_sin;
+ element->E(2*i + 1, 2*i) = tmp_sin;
+ element->E(2*i + 1, 2*i + 1) = tmp_cos;
+ }
+
+ //update C matrix
+ element->C = element->P*element->E.transpose();
+
+ //update fl vector, only if Kl has been initialized
+ if(computeFl)
+ {
+ element->fl = element->Kl*element->pl;
}
}
-//Get physical properties of domain
-void getPhysicalProperties(const std::vector<std::string> & keys, double & E, double & nu, double & rho,
- double & bx, double & by)
+// gets physical properties of FEM_domain: Young's modulus "E", Poisson's ratio "nu",
+// density "rho", volumic force along x "bx", volumic force along y "by".
+void getPhysicalProperties(double & E, double & nu, double & rho, double & bx, double & by)
{
+ std::vector<std::string> keys;
+ gmsh::onelab::getNames(keys, "(Volumic Forces|Materials).+");
for (auto &key : keys)
{
// get corresponding value
@@ -113,6 +263,28 @@ void getPhysicalProperties(const std::vector<std::string> & keys, double & E, do
}
}
+// computes the Hooke's matrix in plane stress "H_matrix" based on the Young's modulus "E" and Poisson's ratio "nu".
+Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu)
+{
+ Eigen::Matrix<double, 3, 3> H_matrix;
+ H_matrix(0,0) = 1.;
+ H_matrix(1,0) = nu;
+ H_matrix(2,0) = 0.;
+ H_matrix(0,1) = nu;
+ H_matrix(1,1) = 1.;
+ H_matrix(2,1) = 0.;
+ H_matrix(0,2) = 0.;
+ H_matrix(1,2) = 0.;
+ H_matrix(2,2) = (1-nu)/2;
+ H_matrix = E/(1-nu*nu)*H_matrix;
+ return H_matrix;
+}
+
+// looping over every element in the FEM_domain (dimension "dim" and entity tags "tags"), filling "full_f_vector",
+// with the local volumic forces computed based on the volumic force ("bx","by") and density "rho" with the help of
+// "nodeIndexMap" mapping each nodeTag to its global index.
+// filling the local stiffness matrices stored in the elementData structure contained in "FEM_kinematicsMap".
+// numerical integration performed using Gauss integration based on "integration_rule".
void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by,
const std::string & integration_rule,
int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap,
@@ -135,7 +307,7 @@ void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double
gmsh::model::mesh::getElements(elementTypes, elementTags,
nodeTags, dim, tags[k]);
- // looping over element types, enlever les déclarations des boucles pr aller plus vite (cfr Boman)
+ // looping over element types
for (std::size_t i = 0; i < elementTypes.size(); ++i)
{
// getting element properties
@@ -169,16 +341,16 @@ void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double
/*---------computing the K matrix (and f vector, only volumic part for now)------*/
// looping over the elements (computing elemental K matrices)
- //double test1 = omp_get_wtime();
- //Eigen::setNbThreads(1);
#ifdef _MSC_VER
// https://github.com/tesseract-ocr/tesseract/issues/3285
// OpenMP 2.0 standard is enforced - we may use /openmp:llvm and keep size_t
#pragma omp parallel for
- for (size_t el = 0; el < elementTags[i].size(); el++){
+ for (size_t el = 0; el < elementTags[i].size(); el++)
+ {
#else
#pragma omp parallel for
- for (std::size_t el = 0; el < elementTags[i].size(); el++){
+ for (std::size_t el = 0; el < elementTags[i].size(); el++)
+ {
#endif
// elemental K matrix and f vector initialization
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> K_matrix(2*numNodes,2*numNodes);
@@ -186,12 +358,13 @@ void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double
K_matrix.setZero();
for (int j = 0; j < 2*numNodes; ++j)f_vector[j] = 0.;
- // looping over the Gauss points (formula of slide 17)
- for (std::size_t j = 0; j < weights.size(); ++j){
+ // looping over the Gauss points
+ for (std::size_t gaussIndex = 0; gaussIndex < weights.size(); ++gaussIndex)
+ {
// converting jacobian to Eigen format
Eigen::Matrix<double, 2, 2> jacob2D;
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
- jacob(&jacobians[9*j+9*weights.size()*el], 3, 3);
+ jacob(&jacobians[9*gaussIndex+9*weights.size()*el], 3, 3);
// 3D to 2D (plane stress case)
jacob2D(0,0) = jacob(0,0);
jacob2D(1,0) = jacob(1,0);
@@ -207,27 +380,47 @@ void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double
// computing the B matrix (formula slide 19) and the f vector (also using gauss integration)
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
-
- computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
+ computeBmatrix(gaussIndex, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
- for (int l = 0; l < numNodes; ++l){ // looping over the number of shape functions
- f_vector[2*l] = f_vector[2*l] + rho*basisFunctions[l+numNodes*j]*bx*determinants[j+weights.size()*el]*weights[j]; //pas oublier le dét !
- f_vector[2*l+1] = f_vector[2*l+1] + rho*basisFunctions[l+numNodes*j]*by*determinants[j+weights.size()*el]*weights[j]; // pas oublier le weights[j] non plus !
+ for (int l = 0; l < numNodes; ++l)
+ { // looping over the number of shape functions
+ f_vector[2*l] = f_vector[2*l] + rho*basisFunctions[l+numNodes*gaussIndex]*bx*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
+ f_vector[2*l+1] = f_vector[2*l+1] + rho*basisFunctions[l+numNodes*gaussIndex]*by*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
}
- K_matrix = K_matrix + B_matrix.transpose()*H_matrix*B_matrix*determinants[j+weights.size()*el]*weights[j];
+ K_matrix = K_matrix + B_matrix.transpose()*H_matrix*B_matrix*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
}
/*------------- ASSEMBLY PROCESS ------------*/
FEM_kinematicsMap[elementTags[i][el]].Kl = K_matrix;
- assembleFlocal(i, el, nodeTags, numNodes, nodeIndexMap, full_f_vector, f_vector);
+ assembleFlocal(el, nodeTags[i], numNodes, nodeIndexMap, full_f_vector, f_vector);
}
}
}
}
-//Assembles the local "f vector" into the full f vector
-void assembleFlocal(const int &i, const int &el,
- const std::vector<std::vector<std::size_t>> & nodeTags,
+// computing the strain displacement "B_matrix" at the "gaussIndex"-th gaussian point of the element
+// consisting of "numNodes" nodes, based on the "basisFunctionsGradient", "determinants" and the inverse of the jacobian
+// matrix "jacobinvtrans".
+void computeBmatrix(const int &gaussIndex, const int &numNodes, const std::vector<double> & basisFunctionsGradient,
+ const std::vector<double> & determinants,
+ const Eigen::Matrix<double, 2, 2> & jacobinvtrans,
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & B_matrix)
+{
+ for (int l = 0; l < numNodes; ++l){ // looping over the number of shape functions
+ B_matrix(0,2*l) = jacobinvtrans(0,0)*basisFunctionsGradient[3*l+numNodes*3*gaussIndex]+jacobinvtrans(0,1)*basisFunctionsGradient[3*l+numNodes*3*gaussIndex+1];
+ B_matrix(1,2*l) = 0.;
+ B_matrix(2,2*l) = jacobinvtrans(1,0)*basisFunctionsGradient[3*l+numNodes*3*gaussIndex]+jacobinvtrans(1,1)*basisFunctionsGradient[3*l+numNodes*3*gaussIndex+1];
+ B_matrix(0,2*l+1) = 0.;
+ B_matrix(1,2*l+1) = B_matrix(2,2*l);
+ B_matrix(2,2*l+1) = B_matrix(0,2*l);
+ }
+}
+
+// assembles the local "f_vector" into the full f vector, "elIndex" is the element index with respect to "nodeTags"
+// having "numNodes" nodes, "nodeTags" gathering all the nodeTags of the current element type in the current entity.
+// "nodeIndexMap" associates a global index to each nodetag.
+void assembleFlocal(const int &elIndex,
+ const std::vector<std::size_t> & nodeTags,
const int &numNodes, std::map<int, int> & nodeIndexMap,
std::vector<double> & full_f_vector,
std::vector<double> & f_vector)
@@ -238,48 +431,19 @@ void assembleFlocal(const int &i, const int &el,
for (int j = 0; j < 2*numNodes; ++j)
{
first_node_index = (int) j/2;
- first_node = nodeTags[i][numNodes*el+first_node_index];
+ first_node = nodeTags[numNodes*elIndex+first_node_index];
my_row = nodeIndexMap[first_node] + j % 2;
#pragma omp atomic update
full_f_vector[my_row] += f_vector[j];
}
}
-/*-----------computing the f vector components related to this particular Neumann BC---------*/
-void computeNeumannBC(const int &numNodes, const int &elTyp, const std::vector<std::vector<std::size_t>> & elementTags,
- const double &tx, const double &ty, std::map<int, int> & nodeIndexMap,
- const std::vector<double> & weights, const std::vector<double> & basisFunctions,
- const std::vector<double> & determinants, const std::vector<std::vector<std::size_t>> & nodeTags,
- std::vector<double> & full_f_vector)
-{
- #pragma omp parallel for
- for (std::size_t el = 0; el < elementTags[elTyp].size(); el++){
- std::vector<double> f_vector(2*numNodes);
- //std::cout << "i am thread: " << omp_get_thread_num() << "\n";
- for (int j = 0; j < 2*numNodes; ++j){
- f_vector[j] = 0.;
- }
- // looping over the Gauss points (formula of slide 17)
- for (std::size_t j = 0; j < weights.size(); ++j){
- // looping over the number of shape functions
- for (int l = 0; l < numNodes; ++l){
- f_vector[2*l] = f_vector[2*l] + basisFunctions[l+numNodes*j]*tx*determinants[j+weights.size()*el]*weights[j]; //pas oublier le dét !
- f_vector[2*l+1] = f_vector[2*l+1] + basisFunctions[l+numNodes*j]*ty*determinants[j+weights.size()*el]*weights[j]; // le weight non plus !
- }
- }
- /*------------- ASSEMBLY PROCESS ------------*/
- for (int j = 0; j < 2*numNodes; ++j){
- int first_node_index = (int) j/2;
- int first_node = nodeTags[elTyp][numNodes*el+first_node_index];
- int my_row = nodeIndexMap[first_node] + j % 2;
- #pragma omp atomic update
- full_f_vector[my_row] += f_vector[j];
- }
- }
-}
-
-/*-----------focus on neumann BC's (surface traction) into full_f_vector------------*/
-void includeNeumannBC(const std::vector<std::string> & keys, const std::string & integration_rule,
+// fills the "full_f_vector" with the surfacic contribution of the neumann boundary conditions based on the
+// "nodeIndexMap" associating a global index to each node tag, the local integration of the surfacic forces is
+// performed using Gauss integration based on "integration rule".
+// "BCkeys" gathers the keys associated to boundary conditions in the .geo file (using SetNumber(...)).
+// "groups" associates a dimension and a tag to each physical group.
+void includeNeumannBC(const std::vector<std::string> & BCkeys, const std::string & integration_rule,
std::map<int, int> & nodeIndexMap,
std::map<std::string, std::pair<int, int>> & groups,
std::vector<double> & full_f_vector)
@@ -300,7 +464,7 @@ void includeNeumannBC(const std::vector<std::string> & keys, const std::string &
int numComponents, numOrientations;
std::vector<double> basisFunctions;
- for (auto &key : keys)
+ for (auto &key : BCkeys)
{
// get corresponding value
std::vector<double> value;
@@ -323,9 +487,8 @@ void includeNeumannBC(const std::vector<std::string> & keys, const std::string &
if(tx || ty) // there is no need the consider Neumann BC if it is homogeneous
groupname = words[1];
}
- //here, if we are dealing with neumann BC's, we should have retrieved a non-empty group_name
+ //here, if we are dealing with neumann BC's, we should have retrieved a non-empty groupname
if(groupname.compare("") != 0){
- //loop over elements of group_name
dim = groups[groupname].first;
tag = groups[groupname].second;
if (tag == 0)
@@ -333,7 +496,7 @@ void includeNeumannBC(const std::vector<std::string> & keys, const std::string &
std::cerr << "Group '" << groupname << "' does not exist!\n";
return;
}
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags); // entities of particular groupname
// loop over entities
for (std::size_t k = 0; k < tags.size(); ++k)
@@ -351,25 +514,44 @@ void includeNeumannBC(const std::vector<std::string> & keys, const std::string &
// get Gauss points coordinates and weights
gmsh::model::mesh::getIntegrationPoints(elementTypes[i],
- integration_rule, //maybe use different integration method on boundary than in domain
+ integration_rule,
localCoords, weights);
// get determinants and Jacobians (only determinants will be useful)
std::vector<double> jacobians, determinants, coords;
gmsh::model::mesh::getJacobians(elementTypes[i],
localCoords, jacobians,
- determinants, coords, tags[k]); // ai dû rajouter mon tags[k] pr focus on this particular entity
+ determinants, coords, tags[k]);
gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords,
"Lagrange", numComponents,
basisFunctions,
numOrientations);
- /*-----------computing the f vector components related to this Neumann BC---------*/
//looping over the element types
- computeNeumannBC(numNodes, i, elementTags, tx, ty, nodeIndexMap,
- weights, basisFunctions,
- determinants, nodeTags, full_f_vector);
+ #pragma omp parallel for
+ for (std::size_t el = 0; el < elementTags[i].size(); el++){
+ std::vector<double> f_vector(2*numNodes);
+ for (int j = 0; j < 2*numNodes; ++j){
+ f_vector[j] = 0.;
+ }
+ // looping over the Gauss points
+ for (std::size_t j = 0; j < weights.size(); ++j){
+ // looping over the number of shape functions
+ for (int l = 0; l < numNodes; ++l){
+ f_vector[2*l] = f_vector[2*l] + basisFunctions[l+numNodes*j]*tx*determinants[j+weights.size()*el]*weights[j]; //pas oublier le dét !
+ f_vector[2*l+1] = f_vector[2*l+1] + basisFunctions[l+numNodes*j]*ty*determinants[j+weights.size()*el]*weights[j]; // le weight non plus !
+ }
+ }
+ /*------------- ASSEMBLY PROCESS ------------*/
+ for (int j = 0; j < 2*numNodes; ++j){
+ int first_node_index = (int) j/2;
+ int first_node = nodeTags[i][numNodes*el+first_node_index];
+ int my_row = nodeIndexMap[first_node] + j % 2;
+ #pragma omp atomic update
+ full_f_vector[my_row] += f_vector[j];
+ }
+ }
}
}
}
@@ -377,364 +559,34 @@ void includeNeumannBC(const std::vector<std::string> & keys, const std::string &
}
}
-// looping over the boundary conditions imposed in the .geo file, the purpose of the whole loop is
-// to fill the two vectors defined just above
-void fillDOFindicesDirBC(const std::vector<std::string> & keys,
- std::map<int, int> & nodeIndexMap,
- std::map<std::string, std::pair<int, int>> & groups,
- std::vector<int> & new_DOFindices,
- std::vector<double> & dir_BC)
+// for the coupled solver, fills the "full_f_vector" with the electrostatic pressure contained in the
+// "electrostaticPressure" map for different element tags, using Gauss integration "integration_rule", the "nodeIndexMap"
+// associating a global nodal index to each nodeTag, "nodeCoordsMap" gathering the coordinates of each node and "groups"
+// containing the different physical groups of the domain.
+void applyElecPressure(const std::string & integration_rule, std::map<int,std::pair<double,double>> & nodeCoordsMap,
+ std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups,
+ std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector)
{
-
- for (std::size_t j = 0; j < new_DOFindices.size(); ++j){
- new_DOFindices[j] = (int) j;
- }
- for (std::size_t j = 0; j < dir_BC.size(); ++j){
- dir_BC[j] = 0;
- }
-
- double ux; // imposed displacement along x axis
- double uy; // imposed displacement along y axis
-
- int dim;
- int tag;
+ // electrostatic pressure must only be applied on the BEM_FEM_boundary
+ std::string groupname = "BEM_FEM_boundary";
+ int dim = groups[groupname].first;
+ int tag = groups[groupname].second;
std::vector<int> tags;
- std::vector<int> elementTypes;
- std::vector<std::vector<std::size_t>> elementTags;
- std::vector<std::vector<std::size_t>> nodeTags;
- std::string elementName;
- int elDim, order, numNodes, numPrimaryNodes;
- std::vector<double> localNodeCoord;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+ int first_nodeTag;
+ int second_nodeTag;
+ double nx;
+ double ny;
+ double p;
+ double tx;
+ double ty;
+ double norm;
int first_node_index;
int first_node;
int my_row;
- for (auto &key : keys)
- {
- // get corresponding value
- std::vector<double> value;
- gmsh::onelab::getNumber(key, value);
- // expected key structure is "type/group_name/field"
- // => split the key string into components
- std::stringstream ss(key);
- std::vector<std::string> words;
- std::string word;
- while(std::getline(ss, word, '/')) // read string until '/'
- words.push_back(word);
- if(words.size()==3)
- {
- if(words[2].compare("ux") == 0 || words[2].compare("uy") == 0){
- ux = 0;
- uy = 0;
- if(words[2].compare("ux") == 0)
- ux = value[0];
- if(words[2].compare("uy") == 0)
- uy = value[0];
- std::string groupname = words[1];
-
- //loop over elements of group_name
- dim = groups[groupname].first;
- tag = groups[groupname].second;
- if (tag == 0)
- {
- std::cerr << "Group '" << groupname << "' does not exist!\n";
- return;
- }
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
-
- for (std::size_t k = 0; k < tags.size(); ++k)
- {
- gmsh::model::mesh::getElements(elementTypes, elementTags,
- nodeTags, dim, tags[k]);
-
- for (std::size_t i = 0; i < elementTypes.size(); ++i)
- {
- gmsh::model::mesh::getElementProperties(elementTypes[i],
- elementName, elDim, order,
- numNodes, localNodeCoord,
- numPrimaryNodes);
- //looping over the elements
- for (std::size_t el = 0; el < elementTags[i].size(); el++){
- for (int j = 0; j < 2*numNodes; ++j){
- first_node_index = (int) j/2;
- first_node = nodeTags[i][numNodes*el+first_node_index];
- my_row = nodeIndexMap[first_node] + j % 2; // j % 2 == 1 corresponds to y nodal displacement
- if(words[2].compare("uy") == 0 && j % 2){ // we are dealing with a uy type of node
- dir_BC[my_row] = uy;
- new_DOFindices[my_row] = -1;
- }
- if(words[2].compare("ux") == 0 && j % 2 == 0){ // ux type of node
- dir_BC[my_row] = ux;
- new_DOFindices[my_row] = -1;
- }
- }
- }
- }
- }
- }
- }
- }
-}
-
-void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups,
- std::map<int, int> & nodeIndexMap, std::vector<double> & final_nodal_forces,
- const std::vector<std::string> & keys)
-{
- double Rx;
- double Ry;
-
- int dim;
- int tag;
- std::vector<int> tags;
-
- std::vector<std::size_t> nodeTags;
- std::vector<double> coord; //useless but mandatory
- std::vector<double> parametricCoord; //useless but mandatory
-
- for (auto &key : keys)
- {
- // get corresponding value
- std::vector<double> value;
- gmsh::onelab::getNumber(key, value);
- // expected key structure is "type/group_name/field"
- // => split the key string into components
- std::stringstream ss(key);
- std::vector<std::string> words;
- std::string word;
- while(std::getline(ss, word, '/')) // read string until '/'
- words.push_back(word);
- if(words.size()==3)
- {
- if(words[2].compare("ux") == 0 || words[2].compare("uy") == 0){
- Rx = 0;
- Ry = 0;
- std::string groupname = words[1];
- std::cout << "-------------------- " << groupname << " ------------------- \n";
-
- //loop over elements of group_name
- dim = groups[groupname].first;
- tag = groups[groupname].second;
- if (tag == 0)
- {
- std::cerr << "Group '" << groupname << "' does not exist!\n";
- return;
- }
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
-
- // loop over entities
- for (std::size_t k = 0; k < tags.size(); ++k)
- {
- // retrieving all the nodes in this particular entity
- gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord, dim, tags[k], true);
-
-
- for (std::size_t j = 0; j < nodeTags.size(); ++j){
- if(words[2].compare("ux") == 0){ // computing Rx
- int my_new_col = nodeIndexMap[nodeTags[j]]; // le noeud numéro nodeIndexMap[nodeTags[j]] correspond à la ligne/colonne dans K, f
- Rx += final_nodal_forces[my_new_col];
- }
- if(words[2].compare("uy") == 0){ // computing Ry
- int my_new_col = nodeIndexMap[nodeTags[j]] + 1;
- Ry += final_nodal_forces[my_new_col];
- }
- }
- }
- std::cout << "Rx: " << Rx << " [N/m], Ry: " << Ry << " [N/m] \n";
- }
- }
- }
-}
-
-void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector, std::vector<std::string> & names)
-{
- int viewtagUx = gmsh::view::add("ux");
- int viewtagUy = gmsh::view::add("uy");
- std::vector<std::vector<double>> dataUx(FEM_nodeTags.size());
- std::vector<std::vector<double>> dataUy(FEM_nodeTags.size());
- for (std::size_t i = 0; i < FEM_nodeTags.size(); ++i)
- {
- dataUx[i].resize(1);
- dataUx[i][0] = full_d_vector[2*i]; //ux // Ã modifier ici
- dataUy[i].resize(1);
- dataUy[i][0] = full_d_vector[2*i+1]; //uy
- }
-
- gmsh::view::addModelData(viewtagUx, 0, names[0], "NodeData",
- FEM_nodeTags, dataUx); //independent of time here
- gmsh::view::addModelData(viewtagUy, 0, names[0], "NodeData",
- FEM_nodeTags, dataUy); //independent of time here
-}
-
-void computeStressStrainElNodal(int & dim, std::vector<int> & tags,
- std::vector<std::size_t> & elems2D, const Eigen::Matrix<double, 3, 3> &H_matrix,
- const double nu, const double E, std::map<int, elementData> & FEM_kinematicsMap,
- std::vector<std::vector<std::vector<double>>> & dataEps, std::vector<std::vector<std::vector<double>>> & dataSig,
- double & max_VM_stress, int & max_VM_nodeTag)
-{
- int data_index = 0;
-
- std::vector<int> elementTypes;
- std::vector<std::vector<std::size_t>> elementTags;
- std::vector<std::vector<std::size_t>> elnodeTags;
- for(std::size_t entityTag = 0; entityTag < tags.size(); ++entityTag)
- {
- gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[entityTag]);
-
- std::string name;
- int dim, order, numNodes, numPrimaryNodes;
- std::vector<double> paramCoord2D;
- int numComponents, numOrientations;
- std::vector<double> basisFunctionsGradient;
-
- Eigen::Vector3d epsilon; // will contain local values of the strains (Voigt notation)
- Eigen::Vector3d sigma; // will contain local values of the stresses (Voigt notation)
-
- Eigen::Matrix<double, 2, 2> jacob2D;
- Eigen::Matrix<double, 2, 2> jacobinv;
- Eigen::Matrix<double, 2, 2> jacobinvtrans;
-
- elementData* element;
- double theta;
-
- Eigen::Matrix<double, 2, 2> strain_tensor;
- Eigen::Matrix<double, 2, 2> stress_tensor;
- Eigen::Matrix<double, 2, 2> rotation_tensor;
-
- for (std::size_t ityp = 0; ityp < elementTypes.size(); ++ityp)
- {
- // get info about this type of element
- gmsh::model::mesh::getElementProperties(elementTypes[ityp], name, dim, order,
- numNodes, paramCoord2D, numPrimaryNodes);
- //paramCoord contains the coordinates of the nodes in 2D
- // passage en 3D pour utiliser getJacobians
- std::vector<double> paramCoord(3*paramCoord2D.size()/2);
- for(std::size_t j = 0; j < paramCoord2D.size()/2; ++j){
- paramCoord[3*j] = paramCoord2D[2*j];
- paramCoord[3*j+1] = paramCoord2D[2*j+1];
- paramCoord[3*j+2] = 0.;
- }
-
- // select element/node tags for this type of element
- auto &etags = elementTags[ityp];
-
- // get determinants and Jacobians (only jacob useful here) AT NODES AND NOT AT GAUSS POINTS
- std::vector<double> jacobians, determinants, coords;
- gmsh::model::mesh::getJacobians(elementTypes[ityp],
- paramCoord, jacobians,
- determinants, coords, tags[entityTag]);
-
- // derivatives of the shape functions at the nodes of the reference element
- // we are here working with derivatives in the reference space, same for all elements
- gmsh::model::mesh::getBasisFunctions(elementTypes[ityp], paramCoord,
- "GradLagrange", numComponents,
- basisFunctionsGradient,
- numOrientations);
-
- // loop over elements of type "ityp"
- for (std::size_t i = 0; i < etags.size(); ++i)
- {
- element = & FEM_kinematicsMap[etags[i]];
- elems2D[data_index] = etags[i];
-
- theta = element->theta;
- rotation_tensor(0,0) = cos(theta);
- rotation_tensor(0,1) = sin(theta);
- rotation_tensor(1,0) = -sin(theta);
- rotation_tensor(1,1) = cos(theta);
-
- dataEps[0][data_index].resize(numNodes);
- dataEps[1][data_index].resize(numNodes);
- dataEps[2][data_index].resize(numNodes);
- dataEps[3][data_index].resize(numNodes);
- dataSig[0][data_index].resize(numNodes);
- dataSig[1][data_index].resize(numNodes);
- dataSig[2][data_index].resize(numNodes);
- dataSig[3][data_index].resize(numNodes);
-
- for (int j = 0; j < numNodes; ++j){
- // converting jacobian to Eigen format
- Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
- jacob(&jacobians[9*j+9*numNodes*i], 3, 3);
- // 3D to 2D (plane stress case)
- jacob2D(0,0) = jacob(0,0);
- jacob2D(1,0) = jacob(1,0);
- jacob2D(0,1) = jacob(0,1);
- jacob2D(1,1) = jacob(1,1);
-
- // computing the inverse
- jacobinv = jacob2D.inverse();
-
- // computing the transpose of the inverse
- jacobinvtrans = jacobinv.transpose();
-
- // computing the B matrix (formula slide 19)
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
-
- computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
-
- // expressed in local axes (and in Voigt notation)
- epsilon = B_matrix*element->pl;
- sigma = H_matrix*epsilon;
-
- // converting to tensor notation (still in local axes)
- strain_tensor(0,0) = epsilon(0);
- strain_tensor(1,0) = epsilon(2)/2; // divide by 2 because engineering shear strain -> pure shear strain
- strain_tensor(0,1) = epsilon(2)/2;
- strain_tensor(1,1) = epsilon(1);
-
- stress_tensor(0,0) = sigma(0);
- stress_tensor(1,0) = sigma(2);
- stress_tensor(0,1) = sigma(2);
- stress_tensor(1,1) = sigma(1);
-
- //converting back to global axes // SUREMENT LE CHANGER !!!!
- strain_tensor = rotation_tensor.transpose()*strain_tensor*rotation_tensor;
- stress_tensor = rotation_tensor.transpose()*stress_tensor*rotation_tensor;
-
- dataEps[0][data_index][j] = strain_tensor(0,0);
- dataEps[1][data_index][j] = strain_tensor(1,1);
- dataEps[2][data_index][j] = strain_tensor(1,0);
- dataEps[3][data_index][j] = -nu/E*(stress_tensor(0,0)+stress_tensor(1,1));// formula slide 11/20 elasticity
- dataSig[0][data_index][j] = stress_tensor(0,0);
- dataSig[1][data_index][j] = stress_tensor(1,1);
- dataSig[2][data_index][j] = stress_tensor(1,0);
- dataSig[3][data_index][j] = sqrt(stress_tensor(0,0)*stress_tensor(0,0) + stress_tensor(1,1)*stress_tensor(1,1) - stress_tensor(0,0)*stress_tensor(1,1) + 3*stress_tensor(1,0)*stress_tensor(1,0));
-
- //retrieving maximal nodal von mises stress
- if(dataSig[3][data_index][j] > max_VM_stress)
- {
- max_VM_stress = dataSig[3][data_index][j];
- max_VM_nodeTag = (int) elnodeTags[ityp][numNodes*i+j];
- }
- }
- data_index++;
- }
- }
- }
-}
-
-void applyElecPressure(const std::string & integration_rule, std::map<int,std::pair<double,double>> & nodeCoordsMap, std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups, std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector)
-{
- std::string groupname = "BEM_FEM_boundary";
- int dim = groups[groupname].first;
- int tag = groups[groupname].second;
- std::vector<int> tags;
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
-
- int first_nodeTag;
- int second_nodeTag;
- double nx;
- double ny;
- double p;
- double tx;
- double ty;
- double norm;
- int first_node_index;
- int first_node;
- int my_row;
-
- for(std::size_t i=0; i< tags.size(); ++i)
+ for(std::size_t i=0; i< tags.size(); ++i)
{
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
@@ -755,14 +607,14 @@ void applyElecPressure(const std::string & integration_rule, std::map<int,std::p
// get Gauss points coordinates and weights
std::vector<double> localCoords, weights;
gmsh::model::mesh::getIntegrationPoints(elementTypes[j],
- integration_rule, //maybe use different integration method on boundary than in domain
+ integration_rule,
localCoords, weights);
// get determinants and Jacobians (only determinants will be useful)
std::vector<double> jacobians, determinants, coords;
gmsh::model::mesh::getJacobians(elementTypes[j],
localCoords, jacobians,
- determinants, coords,tags[i]); // ai dû rajouter mon tags[k] pr focus on this particular entity
+ determinants, coords,tags[i]);
int numComponents, numOrientations;
std::vector<double> basisFunctions;
@@ -775,6 +627,8 @@ void applyElecPressure(const std::string & integration_rule, std::map<int,std::p
for(std::size_t k=0; k < elementTags[j].size(); ++k)
{
std::vector<double> f_vector(2*numNodes);
+ // computing the outward normal to the element, assuming the BEM domain is traveled "aire à gauche", hence
+ // the FEM domain "aire à droite"
first_nodeTag = elnodeTags[j][numNodes*k];
second_nodeTag = elnodeTags[j][numNodes*k+1];
nx = -(nodeCoordsMap[second_nodeTag].second - nodeCoordsMap[first_nodeTag].second);
@@ -783,18 +637,18 @@ void applyElecPressure(const std::string & integration_rule, std::map<int,std::p
nx = nx / norm; //normalisation
ny = ny / norm;
p = electrostaticPressure[elementTags[j][k]];
- tx = p*nx;
+ tx = p*nx; // projecting the pressure onto the outward normal
ty = p*ny;
for (int n = 0; n < 2*numNodes; ++n){
f_vector[n] = 0.;
}
- // looping over the Gauss points (formula of slide 17)
+ // looping over the Gauss points
for (std::size_t n = 0; n < weights.size(); ++n){
// looping over the number of shape functions
for (int l = 0; l < numNodes; ++l){
- f_vector[2*l] = f_vector[2*l] + basisFunctions[l+numNodes*n]*tx*determinants[n+weights.size()*k]*weights[n]; //pas oublier le dét !
- f_vector[2*l+1] = f_vector[2*l+1] + basisFunctions[l+numNodes*n]*ty*determinants[n+weights.size()*k]*weights[n]; // le weight non plus !
+ f_vector[2*l] = f_vector[2*l] + basisFunctions[l+numNodes*n]*tx*determinants[n+weights.size()*k]*weights[n];
+ f_vector[2*l+1] = f_vector[2*l+1] + basisFunctions[l+numNodes*n]*ty*determinants[n+weights.size()*k]*weights[n];
}
}
/*------------- ASSEMBLY PROCESS ------------*/
@@ -810,93 +664,123 @@ void applyElecPressure(const std::string & integration_rule, std::map<int,std::p
}
}
-int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, elementData> & FEM_kinematicsMap, std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags)
+// looping over the boundary conditions imposed in the .geo file and gathered in "BCkeys" related to the
+// physical groups "groups", the goal is to fill the two vectors "new_DOFindices" and "dir_BC":
+// "new_DOFindices" vector will contain the new indices of the K matrix after the "number_removed_lines" rows/columns
+// corresponding to dirichlet boundary conditions have been removed. A removed index will correspond to "new_DOFindices"= -1.
+// dir_BC[index] will contain the dirichlet boundary condition prescribed to the node whose global index is retrieved using
+// "nodeIndexMap".
+void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys,
+ std::map<int, int> & nodeIndexMap,
+ std::map<std::string, std::pair<int, int>> & groups,
+ std::vector<int> & new_DOFindices,
+ std::vector<double> & dir_BC, int & number_removed_lines)
{
+
+ for (std::size_t j = 0; j < new_DOFindices.size(); ++j){
+ new_DOFindices[j] = (int) j;
+ }
+ for (std::size_t j = 0; j < dir_BC.size(); ++j){
+ dir_BC[j] = 0;
+ }
+
+ double ux; // imposed displacement along x axis
+ double uy; // imposed displacement along y axis
+
+ int dim;
+ int tag;
+ std::vector<int> tags;
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
- std::vector<std::vector<std::size_t>> elnodeTags;
-
+ std::vector<std::vector<std::size_t>> nodeTags;
std::string elementName;
int elDim, order, numNodes, numPrimaryNodes;
std::vector<double> localNodeCoord;
-
- int nbelems=0; // will store total number of elements in FEM domain
- for(std::size_t i=0; i< tags.size(); ++i)
+
+ int first_node_index;
+ int first_node;
+ int my_row;
+
+ for (auto &key : BCkeys)
{
- gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[i]);
- for(std::size_t j=0; j< elementTags.size(); ++j)
+ // get corresponding value
+ std::vector<double> value;
+ gmsh::onelab::getNumber(key, value);
+ // expected key structure is "type/group_name/field"
+ // => split the key string into components
+ std::stringstream ss(key);
+ std::vector<std::string> words;
+ std::string word;
+ while(std::getline(ss, word, '/')) // read string until '/'
+ words.push_back(word);
+ if(words.size()==3)
{
- nbelems+=elementTags[j].size();
- // filling the FEM_elementTags vector
- FEM_elementTags.insert(FEM_elementTags.end(), elementTags[j].begin(), elementTags[j].end());
-
- // getting element properties, in particular: number of nodes
- gmsh::model::mesh::getElementProperties(elementTypes[j],
- elementName, elDim, order,
- numNodes, localNodeCoord,
- numPrimaryNodes);
- for(std::size_t k = 0; k < elementTags[j].size(); k++)
- {
- elementData element;
- std::vector<size_t> tmp_nodes(numNodes);
- std::vector<double> tmp_x(numNodes);
- std::vector<double> tmp_y(numNodes);
- std::vector<double> tmp_u(numNodes);
- std::vector<double> tmp_v(numNodes);
- std::vector<double> tmp_ui(numNodes);
- std::vector<double> tmp_vi(numNodes);
- double tmp_xc = 0;
- double tmp_yc = 0;
- for(int l = 0; l < numNodes; l++)
+ if(words[2].compare("ux") == 0 || words[2].compare("uy") == 0){
+ ux = 0;
+ uy = 0;
+ if(words[2].compare("ux") == 0)
+ ux = value[0];
+ if(words[2].compare("uy") == 0)
+ uy = value[0];
+ std::string groupname = words[1];
+
+ //loop over elements of group_name
+ dim = groups[groupname].first;
+ tag = groups[groupname].second;
+ if (tag == 0)
{
- tmp_nodes[l] = elnodeTags[j][numNodes*k + l];
- tmp_x[l] = FEM_nodeCoordsMap[tmp_nodes[l]].first;
- tmp_y[l] = FEM_nodeCoordsMap[tmp_nodes[l]].second;
- tmp_u[l] = 0; // initializing the displacement at zero
- tmp_v[l] = 0;
- tmp_ui[l] = 0;
- tmp_vi[l] = 0;
- tmp_xc += tmp_x[l];
- tmp_yc += tmp_y[l];
+ std::cerr << "Group '" << groupname << "' does not exist!\n";
+ return;
}
- tmp_xc = tmp_xc / numNodes;
- tmp_yc = tmp_yc / numNodes;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
- std::vector<double> tmp_xi(numNodes);
- std::vector<double> tmp_yi(numNodes);
- for(int l = 0; l < numNodes; l++)
+ for (std::size_t k = 0; k < tags.size(); ++k)
{
- tmp_xi[l] = tmp_x[l] - tmp_xc;
- tmp_yi[l] = tmp_y[l] - tmp_yc;
+ gmsh::model::mesh::getElements(elementTypes, elementTags,
+ nodeTags, dim, tags[k]);
+
+ for (std::size_t i = 0; i < elementTypes.size(); ++i)
+ {
+ gmsh::model::mesh::getElementProperties(elementTypes[i],
+ elementName, elDim, order,
+ numNodes, localNodeCoord,
+ numPrimaryNodes);
+ //looping over the elements
+ for (std::size_t el = 0; el < elementTags[i].size(); el++){
+ for (int j = 0; j < 2*numNodes; ++j){
+ first_node_index = (int) j/2;
+ first_node = nodeTags[i][numNodes*el+first_node_index];
+ my_row = nodeIndexMap[first_node] + j % 2; // j % 2 == 1 corresponds to y nodal displacement
+ if(words[2].compare("uy") == 0 && j % 2){ // we are dealing with a uy type of node
+ dir_BC[my_row] = uy;
+ new_DOFindices[my_row] = -1;
+ }
+ if(words[2].compare("ux") == 0 && j % 2 == 0){ // ux type of node
+ dir_BC[my_row] = ux;
+ new_DOFindices[my_row] = -1;
+ }
+ }
+ }
+ }
}
- element.nodes = tmp_nodes;
- element.x = tmp_x;
- element.y = tmp_y;
- element.xc = tmp_xc;
- element.yc = tmp_yc;
- element.xi = tmp_xi;
- element.yi = tmp_yi;
- element.u = tmp_u;
- element.v = tmp_v;
- element.uc = 0;
- element.vc = 0;
- element.theta = 0;
- element.ui = tmp_ui;
- element.vi = tmp_vi;
- element.pl = Eigen::Matrix<double, Eigen::Dynamic, 1>::Zero(2*numNodes,1);
- element.A = Eigen::Matrix<double, Eigen::Dynamic, 1>(2*numNodes,1);
- element.G = Eigen::Matrix<double, 1, Eigen::Dynamic>(1,2*numNodes);
- element.E = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Zero(2*numNodes, 2*numNodes);
-
- FEM_kinematicsMap[elementTags[j][k]] = element;
- updateKinMatrices(& FEM_kinematicsMap[elementTags[j][k]], false);
}
}
}
- return nbelems;
+ for (std::size_t j = 0; j < new_DOFindices.size(); ++j){
+ if(new_DOFindices[j] < 0){ // line has to be removed
+ number_removed_lines++;
+ }
+ else{
+ new_DOFindices[j] -= number_removed_lines;
+ }
+ }
}
-void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector, std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags)
+// updating the kinematics of the "FEM_elementTags" contained in the "FEM_kinematicsMap",
+// based on the "full_d_vector" gathering the nodal displacements associated to the global
+// node indices contained in "nodeIndexMap".
+void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector,
+ std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags)
{
size_t numNodes;
size_t tmp_nodeTag;
@@ -946,7 +830,7 @@ void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vecto
{
angle2 = angle1 - PI;
}
- // filling ui and vi while finding the angle minimizing equation (4) in the paper (what i call the norm)
+ // filling ui and vi while finding the angle minimizing equation (70) in the report (VERIFIER EQUATION)
std::vector<double> tmp_ui1(numNodes);
std::vector<double> tmp_ui2(numNodes);
std::vector<double> tmp_vi1(numNodes);
@@ -965,76 +849,67 @@ void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vecto
if(norm1 < norm2)
{
element->theta = angle1;
- element->ui = tmp_ui1;
- element->vi = tmp_vi1;
- }
- else
- {
- element->theta = angle2;
- element->ui = tmp_ui2;
- element->vi = tmp_vi2;
- }
- for(size_t l = 0; l < numNodes; l++)
- {
- element->pl(2*l) = element->ui[l];
- element->pl(2*l+1) = element->vi[l];
- }
- updateKinMatrices(element, true);
- }
-}
-
-void updateKinMatrices(elementData* element, bool computeFl)
-{
- size_t numNodes = element->nodes.size();
- // update A matrix
- for(size_t i = 0; i < numNodes; i++)
- {
- element->A(2*i, 0) = -(element->vi[i]+element->yi[i]);
- element->A(2*i + 1, 0) = element->ui[i]+element->xi[i];
- }
-
- // update G matrix
- double denominator = 0;
- for(size_t i = 0; i < numNodes; i++)
- {
- element->G(0, 2*i) = -element->yi[i];
- element->G(0, 2*i + 1) = element->xi[i];
- denominator += element->xi[i]*(element->ui[i]+element->xi[i]) + element->yi[i]*(element->vi[i]+element->yi[i]);
+ element->ui = tmp_ui1;
+ element->vi = tmp_vi1;
+ }
+ else
+ {
+ element->theta = angle2;
+ element->ui = tmp_ui2;
+ element->vi = tmp_vi2;
+ }
+ for(size_t l = 0; l < numNodes; l++)
+ {
+ element->pl(2*l) = element->ui[l];
+ element->pl(2*l+1) = element->vi[l];
+ }
+ updateKinMatrices(element, true);
}
- element->G *= 1 / denominator;
-
- // update P matrix
- element->P = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Identity(2*numNodes,2*numNodes) - element->A*element->G;
+}
- // update E matrix
- double tmp_cos = cos(element->theta);
- double tmp_sin = sin(element->theta);
- for(size_t i = 0; i < numNodes; i++)
+// retrieving the global nodal forces based on local nodal forces of the
+// "FEM_elementTags" contained in the "FEM_kinematicsMap",
+// each node Tag corresponds to a global index given by "nodeIndexMap".
+void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags)
+{
+ for(size_t i = 0; i < final_nodal_forces.size(); i++)
{
- element->E(2*i, 2*i) = tmp_cos;
- element->E(2*i, 2*i + 1) = -tmp_sin;
- element->E(2*i + 1, 2*i) = tmp_sin;
- element->E(2*i + 1, 2*i + 1) = tmp_cos;
+ final_nodal_forces[i] = 0;
}
- //update B matrix
- element->B = element->P*element->E.transpose();
-
- //update fl vector, only if Kl has been initialized
- if(computeFl)
+ #pragma omp parallel for
+ for(size_t i = 0; i < FEM_elementTags.size(); i++)
{
- element->fl = element->Kl*element->pl;
+ elementData* element = & FEM_kinematicsMap[FEM_elementTags[i]];
+ int numNodes = (int) element->nodes.size();
+ Eigen::Matrix<double, Eigen::Dynamic, 1> fg = element->C.transpose()*element->fl;
+ for(int j = 0; j < numNodes; j++)
+ {
+ int nodeTag = element->nodes[j];
+ #pragma omp atomic update
+ final_nodal_forces[nodeIndexMap[nodeTag]] += fg(2*j);
+ #pragma omp atomic update
+ final_nodal_forces[nodeIndexMap[nodeTag] + 1] += fg(2*j + 1);
+ }
}
}
-void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elementData> & FEM_kinematicsMap, std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, std::vector<std::size_t> & FEM_elementTags)
+// assembles the "global_f_vector" based on the local nodal forces vector of the
+// "FEM_elementTags" contained in the "FEM_kinematicsMap",
+// without taking into account the DOFs fixed by a dirichlet boundary condition,
+// the nodes are renumbered using "new_DOFindices".
+// each node Tag corresponds to a global index given by "nodeIndexMap".
+void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices,
+ std::vector<std::size_t> & FEM_elementTags)
{
#pragma omp parallel for
for(size_t i = 0; i < FEM_elementTags.size(); i++)
{
elementData* element = & FEM_kinematicsMap[FEM_elementTags[i]];
int numNodes = (int) element->nodes.size();
- Eigen::Matrix<double, Eigen::Dynamic, 1> fg = element->B.transpose()*element->fl;
+ Eigen::Matrix<double, Eigen::Dynamic, 1> fg = element->C.transpose()*element->fl;
for(int j = 0; j < numNodes; j++)
{
int nodeTag = element->nodes[j];
@@ -1054,21 +929,13 @@ void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elemen
}
}
-double computeResidualNorm(std::vector<double> & global_f_vector, std::vector<double> & f_app_vector, Eigen::Matrix<double, Eigen::Dynamic, 1> & residual)
-{
- // norm computed as a geometrical norm (one could also use an other norm)
- double res = 0;
- for(size_t i = 0; i < global_f_vector.size(); i++)
- {
- residual(i,0) = global_f_vector[i] - f_app_vector[i];
- res += residual(i,0)*residual(i,0);
- }
- res = res / f_app_vector.size();
- res = sqrt(res);
- return res;
-}
-
-void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::map<int, elementData> & FEM_kinematicsMap, std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, std::vector<std::size_t> & FEM_elementTags)
+// assembling the global tangent stiffness matrix based on the local elemental stiffness matrices contained
+// in the "FEM_kinematicsMap" associated to each "FEM_elementTags".
+// "new_DOFindices" maps the global initial node index (from "nodeIndexMap") to the index
+// obtained without considering the fixed DOFs.
+void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices,
+ std::vector<std::size_t> & FEM_elementTags)
{
// to allow efficient parallelization, one triplet list is created for every thread, then the different lists are assembled
// to further increase efficiency, space is reserved for each list
@@ -1100,9 +967,10 @@ void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::
F(0,2*j+1) = -tmp_F(2*j, 0);
}
+ // see report for the formula
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Kh = element->E*(-F.transpose()*element->G - element->G.transpose()*F*element->P)*element->E.transpose();
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Kg = element->B.transpose()*element->Kl*element->B + Kh;
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Kg = element->C.transpose()*element->Kl*element->C + Kh;
int threadNum = omp_get_thread_num();
my_Kg_list = & thread_Kg_list[threadNum];
@@ -1145,7 +1013,10 @@ void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::
global_tangent_matrix.setFromTriplets(Kg_tripletList.begin(), Kg_tripletList.end());
}
-void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<double> & global_p_vector, std::vector<int> & new_DOFindices)
+// updating the "full_d_vector" based on the "global_p_vector", using the global and reduced index mapping:
+// "new_DOFindices".
+void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<double> & global_p_vector,
+ std::vector<int> & new_DOFindices)
{
for(size_t i = 0; i < full_d_vector.size(); i++)
{
@@ -1156,99 +1027,350 @@ void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<doub
}
}
-void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map<int, elementData> & FEM_kinematicsMap, std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags)
-{
- for(size_t i = 0; i < final_nodal_forces.size(); i++)
- {
- final_nodal_forces[i] = 0;
- }
+// retrieving the nodal displacements on the BEM_FEM_boundary in order to communicate them to the
+// BEM solver in the two-way coupled solver. Using the "groups" to find the physical group associated
+// to the boundary, all nodes of the boundary are associated through the "boundaryDisplacementMap" their
+// nodal displacement (u,v) retrieved from "full_d_vector" using the indexation of "nodeIndexMap".
+void retrieveBoundaryDisplacements(std::map<int,std::pair<double,double>> & boundaryDisplacementMap,
+ const std::vector<double> & full_d_vector, std::map<int, int> & nodeIndexMap,
+ std::map<std::string, std::pair<int, int>> & groups)
+{
+ std::string groupname = "BEM_FEM_boundary";
+ int dim = groups[groupname].first;
+ int tag = groups[groupname].second;
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> elnodeTags;
+
+ int nodeTag;
+ double u;
+ double v;
+
+ std::string elementName;
+ int elDim, order, numNodes, numPrimaryNodes;
+ std::vector<double> localNodeCoord;
+
+ for(std::size_t i=0; i< tags.size(); ++i)
+ {
+ gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[i]);
+ for(std::size_t j=0; j< elementTypes.size(); ++j)
+ {
+
+ // getting element properties
+ gmsh::model::mesh::getElementProperties(elementTypes[j],
+ elementName, elDim, order,
+ numNodes, localNodeCoord,
+ numPrimaryNodes);
+
+ for(std::size_t k=0; k < elementTags[j].size(); ++k)
+ {
+ for(int l = 0; l < numNodes; l++)
+ {
+ nodeTag = elnodeTags[j][numNodes*k + l];
+ u = full_d_vector[nodeIndexMap[nodeTag]];
+ v = full_d_vector[nodeIndexMap[nodeTag] + 1];
+ boundaryDisplacementMap[nodeTag] = std::pair<double, double>(u, v);
+ }
+ }
+ }
+ }
+}
+
+// computing the reaction forces on all the physical groups "groups" for which a dirichlet boundary
+// condition has been prescribed (using the "BCkeys" to retrieve the groups from the .geo file),
+// based on the "final_nodal_forces" and the "nodeIndexMap" to retrieve the global index associated
+// to one given node tag.
+void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups,
+ std::map<int, int> & nodeIndexMap, std::vector<double> & final_nodal_forces,
+ const std::vector<std::string> & BCkeys)
+{
+ double Rx;
+ double Ry;
+
+ int dim;
+ int tag;
+ std::vector<int> tags;
+
+ std::vector<std::size_t> nodeTags;
+ std::vector<double> coord; //useless but mandatory
+ std::vector<double> parametricCoord; //useless but mandatory
+
+ for (auto &key : BCkeys)
+ {
+ // get corresponding value
+ std::vector<double> value;
+ gmsh::onelab::getNumber(key, value);
+ // expected key structure is "type/group_name/field"
+ // => split the key string into components
+ std::stringstream ss(key);
+ std::vector<std::string> words;
+ std::string word;
+ while(std::getline(ss, word, '/')) // read string until '/'
+ words.push_back(word);
+ if(words.size()==3)
+ {
+ if(words[2].compare("ux") == 0 || words[2].compare("uy") == 0){
+ Rx = 0;
+ Ry = 0;
+ std::string groupname = words[1];
+ std::cout << "-------------------- " << groupname << " ------------------- \n";
+
+ //loop over elements of group_name
+ dim = groups[groupname].first;
+ tag = groups[groupname].second;
+ if (tag == 0)
+ {
+ std::cerr << "Group '" << groupname << "' does not exist!\n";
+ return;
+ }
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+ // loop over entities
+ for (std::size_t k = 0; k < tags.size(); ++k)
+ {
+ // retrieving all the nodes in this particular entity
+ gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord, dim, tags[k], true);
+
+
+ for (std::size_t j = 0; j < nodeTags.size(); ++j){
+ if(words[2].compare("ux") == 0){ // computing Rx
+ int my_new_col = nodeIndexMap[nodeTags[j]]; //nodeIndexMap[nodeTags[j]] corresponds to line/column index in K, f
+ Rx += final_nodal_forces[my_new_col];
+ }
+ if(words[2].compare("uy") == 0){ // computing Ry
+ int my_new_col = nodeIndexMap[nodeTags[j]] + 1;
+ Ry += final_nodal_forces[my_new_col];
+ }
+ }
+ }
+ std::cout << "Rx: " << Rx << " [N/m], Ry: " << Ry << " [N/m] \n";
+ }
+ }
+ }
+}
+
+// plotting the nodal displacement values associated to "FEM_nodeTags" in the "names[0]" model,
+// based on the "full_d_vector".
+void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector,
+ std::vector<std::string> & names)
+{
+ int viewtagUx = gmsh::view::add("ux");
+ int viewtagUy = gmsh::view::add("uy");
+ std::vector<std::vector<double>> dataUx(FEM_nodeTags.size());
+ std::vector<std::vector<double>> dataUy(FEM_nodeTags.size());
+ for (std::size_t i = 0; i < FEM_nodeTags.size(); ++i)
+ {
+ dataUx[i].resize(1);
+ dataUx[i][0] = full_d_vector[2*i]; //ux // Ã modifier ici
+ dataUy[i].resize(1);
+ dataUy[i][0] = full_d_vector[2*i+1]; //uy
+ }
+
+ gmsh::view::addModelData(viewtagUx, 0, names[0], "NodeData",
+ FEM_nodeTags, dataUx); //independent of time here
+ gmsh::view::addModelData(viewtagUy, 0, names[0], "NodeData",
+ FEM_nodeTags, dataUy); //independent of time here
+}
+
+// computing the elemental-nodal strains and stresses in the FEM domain (dimension "dim" and entity tags "tags"),
+// based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", looping over the
+// "nbelems" elements contained in "FEM_kinematicsMap".
+// returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag".
+void computeStressStrainElNodal(int & dim, std::vector<int> & tags, const Eigen::Matrix<double, 3, 3> &H_matrix,
+ const double nu, const double E, std::map<int, elementData> & FEM_kinematicsMap,
+ int & nbelems, double & max_VM_stress, int & max_VM_nodeTag)
+{
+ std::vector<std::string> names;
+ gmsh::model::list(names);
+
+ int viewTagEpsX = gmsh::view::add("nodal_eps_x");
+ int viewTagEpsY = gmsh::view::add("nodal_eps_y");
+ int viewTagEpsXY = gmsh::view::add("nodal_2eps_xy");
+ int viewTagEpsZ = gmsh::view::add("nodal_eps_z");
+ int viewTagSigX = gmsh::view::add("nodal_sig_x");
+ int viewTagSigY = gmsh::view::add("nodal_sig_y");
+ int viewTagSigXY = gmsh::view::add("nodal_sig_xy");
+ int viewTagSigVM = gmsh::view::add("nodal_sig_VM");
+
+ std::vector<std::vector<std::vector<double>>> dataEps(4);
+ dataEps[0].resize(nbelems); // eps X
+ dataEps[1].resize(nbelems); // eps Y
+ dataEps[2].resize(nbelems); // eps XY
+ dataEps[3].resize(nbelems); // eps Z
+
+ std::vector<std::vector<std::vector<double>>> dataSig(4);
+ dataSig[0].resize(nbelems); // Sig X
+ dataSig[1].resize(nbelems); // Sig Y
+ dataSig[2].resize(nbelems); // Sig XY
+ dataSig[3].resize(nbelems); // Sig VM
+ std::vector<std::size_t> elems2D(nbelems);
+
+ int data_index = 0;
+
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> elnodeTags;
+ for(std::size_t entityTag = 0; entityTag < tags.size(); ++entityTag)
+ {
+ gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[entityTag]);
+
+ std::string name;
+ int dim, order, numNodes, numPrimaryNodes;
+ std::vector<double> paramCoord2D;
+ int numComponents, numOrientations;
+ std::vector<double> basisFunctionsGradient;
+
+ Eigen::Vector3d epsilon; // will contain local values of the strains (Voigt notation)
+ Eigen::Vector3d sigma; // will contain local values of the stresses (Voigt notation)
+
+ Eigen::Matrix<double, 2, 2> jacob2D;
+ Eigen::Matrix<double, 2, 2> jacobinv;
+ Eigen::Matrix<double, 2, 2> jacobinvtrans;
+
+ elementData* element;
+ double theta;
+
+ Eigen::Matrix<double, 2, 2> strain_tensor;
+ Eigen::Matrix<double, 2, 2> stress_tensor;
+ Eigen::Matrix<double, 2, 2> rotation_tensor;
+
+ for (std::size_t ityp = 0; ityp < elementTypes.size(); ++ityp)
+ {
+ // get info about this type of element
+ gmsh::model::mesh::getElementProperties(elementTypes[ityp], name, dim, order,
+ numNodes, paramCoord2D, numPrimaryNodes);
+ //paramCoord contains the coordinates of the nodes in 2D
+ // passage en 3D pour utiliser getJacobians
+ std::vector<double> paramCoord(3*paramCoord2D.size()/2);
+ for(std::size_t j = 0; j < paramCoord2D.size()/2; ++j){
+ paramCoord[3*j] = paramCoord2D[2*j];
+ paramCoord[3*j+1] = paramCoord2D[2*j+1];
+ paramCoord[3*j+2] = 0.;
+ }
+
+ // select element/node tags for this type of element
+ auto &etags = elementTags[ityp];
+
+ // get determinants and Jacobians (only jacob useful here) AT NODES AND NOT AT GAUSS POINTS
+ std::vector<double> jacobians, determinants, coords;
+ gmsh::model::mesh::getJacobians(elementTypes[ityp],
+ paramCoord, jacobians,
+ determinants, coords, tags[entityTag]);
+
+ // derivatives of the shape functions at the nodes of the reference element
+ // we are here working with derivatives in the reference space, same for all elements
+ gmsh::model::mesh::getBasisFunctions(elementTypes[ityp], paramCoord,
+ "GradLagrange", numComponents,
+ basisFunctionsGradient,
+ numOrientations);
+
+ // loop over elements of type "ityp"
+ for (std::size_t i = 0; i < etags.size(); ++i)
+ {
+ element = & FEM_kinematicsMap[etags[i]];
+ elems2D[data_index] = etags[i];
+
+ theta = element->theta;
+ rotation_tensor(0,0) = cos(theta);
+ rotation_tensor(0,1) = sin(theta);
+ rotation_tensor(1,0) = -sin(theta);
+ rotation_tensor(1,1) = cos(theta);
+
+ dataEps[0][data_index].resize(numNodes);
+ dataEps[1][data_index].resize(numNodes);
+ dataEps[2][data_index].resize(numNodes);
+ dataEps[3][data_index].resize(numNodes);
+ dataSig[0][data_index].resize(numNodes);
+ dataSig[1][data_index].resize(numNodes);
+ dataSig[2][data_index].resize(numNodes);
+ dataSig[3][data_index].resize(numNodes);
+
+ for (int j = 0; j < numNodes; ++j){
+ // converting jacobian to Eigen format
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
+ jacob(&jacobians[9*j+9*numNodes*i], 3, 3);
+ // 3D to 2D (plane stress case)
+ jacob2D(0,0) = jacob(0,0);
+ jacob2D(1,0) = jacob(1,0);
+ jacob2D(0,1) = jacob(0,1);
+ jacob2D(1,1) = jacob(1,1);
+
+ // computing the inverse
+ jacobinv = jacob2D.inverse();
- #pragma omp parallel for
- for(size_t i = 0; i < FEM_elementTags.size(); i++)
- {
- elementData* element = & FEM_kinematicsMap[FEM_elementTags[i]];
- int numNodes = (int) element->nodes.size();
- Eigen::Matrix<double, Eigen::Dynamic, 1> fg = element->B.transpose()*element->fl;
- for(int j = 0; j < numNodes; j++)
- {
- int nodeTag = element->nodes[j];
- #pragma omp atomic update
- final_nodal_forces[nodeIndexMap[nodeTag]] += fg(2*j);
- #pragma omp atomic update
- final_nodal_forces[nodeIndexMap[nodeTag] + 1] += fg(2*j + 1);
- }
- }
-}
+ // computing the transpose of the inverse
+ jacobinvtrans = jacobinv.transpose();
-void retrieveBoundaryDisplacements(std::map<int,std::pair<double,double>> & boundaryDisplacementMap, const std::vector<double> & full_d_vector, std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups)
-{
- std::string groupname = "BEM_FEM_boundary";
- int dim = groups[groupname].first;
- int tag = groups[groupname].second;
- std::vector<int> tags;
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
- std::vector<int> elementTypes;
- std::vector<std::vector<std::size_t>> elementTags;
- std::vector<std::vector<std::size_t>> elnodeTags;
+ // computing the B matrix (formula slide 19)
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
- int nodeTag;
- double u;
- double v;
+ computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
+
+ // expressed in local axes (and in Voigt notation)
+ epsilon = B_matrix*element->pl;
+ sigma = H_matrix*epsilon;
- std::string elementName;
- int elDim, order, numNodes, numPrimaryNodes;
- std::vector<double> localNodeCoord;
+ // converting to tensor notation (still in local axes)
+ strain_tensor(0,0) = epsilon(0);
+ strain_tensor(1,0) = epsilon(2)/2; // divide by 2 because engineering shear strain -> pure shear strain
+ strain_tensor(0,1) = epsilon(2)/2;
+ strain_tensor(1,1) = epsilon(1);
- for(std::size_t i=0; i< tags.size(); ++i)
- {
- gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[i]);
- for(std::size_t j=0; j< elementTypes.size(); ++j)
- {
+ stress_tensor(0,0) = sigma(0);
+ stress_tensor(1,0) = sigma(2);
+ stress_tensor(0,1) = sigma(2);
+ stress_tensor(1,1) = sigma(1);
- // getting element properties
- gmsh::model::mesh::getElementProperties(elementTypes[j],
- elementName, elDim, order,
- numNodes, localNodeCoord,
- numPrimaryNodes);
+ strain_tensor = rotation_tensor.transpose()*strain_tensor*rotation_tensor;
+ stress_tensor = rotation_tensor.transpose()*stress_tensor*rotation_tensor;
- for(std::size_t k=0; k < elementTags[j].size(); ++k)
- {
- for(int l = 0; l < numNodes; l++)
- {
- nodeTag = elnodeTags[j][numNodes*k + l];
- u = full_d_vector[nodeIndexMap[nodeTag]];
- v = full_d_vector[nodeIndexMap[nodeTag] + 1];
- boundaryDisplacementMap[nodeTag] = std::pair<double, double>(u, v);
+ dataEps[0][data_index][j] = strain_tensor(0,0);
+ dataEps[1][data_index][j] = strain_tensor(1,1);
+ dataEps[2][data_index][j] = strain_tensor(1,0);
+ dataEps[3][data_index][j] = -nu/E*(stress_tensor(0,0)+stress_tensor(1,1));// from plane stress hypothesis
+ dataSig[0][data_index][j] = stress_tensor(0,0);
+ dataSig[1][data_index][j] = stress_tensor(1,1);
+ dataSig[2][data_index][j] = stress_tensor(1,0);
+ // VM stress in plane stress
+ dataSig[3][data_index][j] = sqrt(stress_tensor(0,0)*stress_tensor(0,0) + stress_tensor(1,1)*stress_tensor(1,1) - stress_tensor(0,0)*stress_tensor(1,1) + 3*stress_tensor(1,0)*stress_tensor(1,0));
+
+ //retrieving maximal nodal von mises stress
+ if(dataSig[3][data_index][j] > max_VM_stress)
+ {
+ max_VM_stress = dataSig[3][data_index][j];
+ max_VM_nodeTag = (int) elnodeTags[ityp][numNodes*i+j];
+ }
}
+ data_index++;
}
}
}
-}
-
-void displayTime(double &start, double &end, std::string &function_name, bool &display_time)
-{
- //displays the time of the desired code section
- if(display_time)
- {
- std::cout << std::endl;
- std::cout << function_name << std::endl;
- std::cout << "Elapsed time in nanoseconds: "
- << end - start << " ns" << std::endl;
-
- std::cout << "Elapsed time in microseconds: "
- << end - start << " µs" << std::endl;
-
- std::cout << "Elapsed time in milliseconds: "
- << end - start << " ms" << std::endl;
- std::cout << "Elapsed time in seconds: "
- << end - start << " sec";
- std::cout << std::endl;
- std::cout << std::endl;
- }
+ gmsh::view::addModelData(viewTagEpsX, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[0], 1);
+ gmsh::view::addModelData(viewTagEpsY, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[1], 1);
+ gmsh::view::addModelData(viewTagEpsXY, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[2], 1);
+ gmsh::view::addModelData(viewTagEpsZ, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[3], 1);
+ gmsh::view::addModelData(viewTagSigX, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[0], 1);
+ gmsh::view::addModelData(viewTagSigY, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[1], 1);
+ gmsh::view::addModelData(viewTagSigXY, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[2], 1);
+ gmsh::view::addModelData(viewTagSigVM, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[3], 1);
}
-void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & full_d_vector, std::vector<double> & final_nodal_forces)
+// computing the work done by external forces, as the nodal forces are energy consistent, it
+// is sufficient to compute the "externalWork" as the dot product of the nodal displacements
+// "full_d_vector" and the nodal forces "final_nodal_forces".
+void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & full_d_vector,
+ std::vector<double> & final_nodal_forces)
{
for(size_t i = 0; i < full_d_vector.size(); i++)
{
@@ -1256,6 +1378,10 @@ void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> &
}
}
+// computing the total (internal) strain energy "strain_energy" in the whole FEM domain (dimension "dim" and
+// entity tags "tags"), using the elementData stored in "FEM_kinematicsMap" gathering the informations for all
+// FEM elements. The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains.
+// Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule".
void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & tags,
std::map<int, elementData> & FEM_kinematicsMap, const Eigen::Matrix<double, 3, 3> &H_matrix,
const std::string & integration_rule)
@@ -1271,9 +1397,6 @@ void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & ta
std::vector<double> basisFunctions;
std::vector<double> basisFunctionsGradient;
- Eigen::Vector3d epsilon; // will contain local values of the strains
- Eigen::Vector3d sigma; // will contain local values of the stresses
-
for (std::size_t k = 0; k < tags.size(); ++k)
{
gmsh::model::mesh::getElements(elementTypes, elementTags,
@@ -1311,19 +1434,19 @@ void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & ta
basisFunctionsGradient,
numOrientations);
-
+ #pragma omp parallel for
for (std::size_t el = 0; el < elementTags[i].size(); el++){
Eigen::Matrix<double, Eigen::Dynamic, 1> local_nodal_displacement = FEM_kinematicsMap[elementTags[i][el]].pl;
double elementalStrainEnergy = 0;
- // looping over the Gauss points (formula of slide 17)
- for (std::size_t j = 0; j < weights.size(); ++j){
+ // looping over the Gauss points: Gauss integration
+ for (std::size_t gaussIndex = 0; gaussIndex < weights.size(); ++gaussIndex){
// converting jacobian to Eigen format
Eigen::Matrix<double, 2, 2> jacob2D;
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
- jacob(&jacobians[9*j+9*weights.size()*el], 3, 3);
+ jacob(&jacobians[9*gaussIndex+9*weights.size()*el], 3, 3);
// 3D to 2D (plane stress case)
jacob2D(0,0) = jacob(0,0);
jacob2D(1,0) = jacob(1,0);
@@ -1336,14 +1459,13 @@ void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & ta
// computing the transpose of the inverse
Eigen::Matrix<double, 2, 2> jacobinvtrans = jacobinv.transpose();
- // computing the B matrix (formula slide 19) and the f vector (also using gauss integration)
+ // computing the B matrix
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
+ computeBmatrix(gaussIndex, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
- computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
-
- epsilon = B_matrix*local_nodal_displacement;
- sigma = H_matrix*epsilon;
+ Eigen::Vector3d epsilon = B_matrix*local_nodal_displacement;
+ Eigen::Vector3d sigma = H_matrix*epsilon;
double dot_product = 0;
for(int z = 0; z < 3; z++)
@@ -1351,67 +1473,58 @@ void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & ta
dot_product += epsilon(z)*sigma(z);
}
- elementalStrainEnergy += 0.5*dot_product*determinants[j+weights.size()*el]*weights[j];
+ elementalStrainEnergy += 0.5*dot_product*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
}
/*------------- ASSEMBLY PROCESS ------------*/
+ #pragma omp atomic update
strainEnergy += elementalStrainEnergy;
}
}
}
}
-/*-------------- FUNCTIONS SPECIFIC TO LINEAR FEM ------------------*/
-
-void assembleKtriplet(const int i, const int el,
- const std::vector<std::vector<std::size_t>> & nodeTags,
- const int numNodes, std::map<int, int> & nodeIndexMap,
- std::vector<Eigen::Triplet <double>> & k_tripletList,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix,
- std::vector<double> & f_vector)
+//displays the time of the desired code section "function_name", if "display_time" is set as true.
+void displayTime(double &start, double &end, std::string &function_name, bool &display_time)
{
- for (int j = 0; j < 2*numNodes; ++j){
- int first_node_index = (int) j/2;
- int first_node = nodeTags[i][numNodes*el+first_node_index];
- int my_row = nodeIndexMap[first_node] + j % 2; // node list begins at 1, '-1' correction to begin at 0
- for(int l = 0; l <= j; ++l){ // to exploit symmetry
- int second_node_index = (int) l/2;
- int second_node = nodeTags[i][numNodes*el+second_node_index];
- int my_col = nodeIndexMap[second_node] + l % 2;
+ if(display_time)
+ {
+ std::cout << std::endl;
+ std::cout << function_name << std::endl;
+ std::cout << "Elapsed time in nanoseconds: "
+ << end - start << " ns" << std::endl;
- k_tripletList.push_back(Eigen::Triplet<double>(my_row,my_col,K_matrix(j,l)));
- if(j != l){
- k_tripletList.push_back(Eigen::Triplet<double>(my_col,my_row,K_matrix(j,l)));
- }
- }
- }
-}
+ std::cout << "Elapsed time in microseconds: "
+ << end - start << " µs" << std::endl;
-//assembles the full matrix "f" of one given thread
-void assembleFthread(const int &i, const int &el,
- const std::vector<std::vector<std::size_t>> & nodeTags,
- const int &numNodes, std::map<int, int> & nodeIndexMap,
- std::vector<std::pair<int,double>> & thread_f_vector,
- std::vector<double> & f_vector)
-{
- int first_node_index;
- int first_node;
- int my_row;
- for (int j = 0; j < 2*numNodes; ++j)
- {
- first_node_index = (int) j/2;
- first_node = nodeTags[i][numNodes*el+first_node_index];
- my_row = nodeIndexMap[first_node] + j % 2;
- thread_f_vector.push_back(std::pair<int, double>(my_row, f_vector[j]));
+ std::cout << "Elapsed time in milliseconds: "
+ << end - start << " ms" << std::endl;
+
+ std::cout << "Elapsed time in seconds: "
+ << end - start << " sec";
+ std::cout << std::endl;
+ std::cout << std::endl;
}
}
+/*-------------- FUNCTIONS SPECIFIC TO LINEAR FEM ------------------*/
-void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by,
- const std::string & integration_rule,
- int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap,
- std::vector<std::vector<Eigen::Triplet <double>>> & my_k_list,
- std::vector<std::vector<std::pair<int,double>>> & thread_f_vector)
+// assembling the full stiffness matrix "full_K_matrix" using first a triplet vector, assembling the volumic part
+// of the nodal forces "full_f_vector", in the FEM_domain (dimension "dim" and entity tag "tags") using the global
+// node indexation defined by "nodeIndexMap". "H_matrix" is the plane stress Hooke's matrix, "rho" the density of the
+// material, "bx" the volumic force along x, "by" the volumic force along y.
+// Numerical integration using gauss integration "integration_rule".
+void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by,
+ const std::string & integration_rule, int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap,
+ Eigen::SparseMatrix<double> & full_K_matrix, std::vector<double> & full_f_vector)
{
+ // the full K matrix will first be assembled as a triplet list (vector) (to increase efficiency),
+ // before initializing the full sparse matrix in one single time;
+ int max_threads = omp_get_max_threads();
+ // one triplet vector declared for every thread (for the K triplet list) (will be gathered at the end of the method)
+ std::vector<std::vector<Eigen::Triplet <double>>> my_k_list(max_threads);
+ // one pair vector declared for every thread (for the f vector) (will be gathered at the end of the method)
+ std::vector<std::vector<std::pair<int,double>>> thread_f_vector(max_threads);
+
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
std::vector<std::vector<std::size_t>> nodeTags;
@@ -1428,7 +1541,7 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
gmsh::model::mesh::getElements(elementTypes, elementTags,
nodeTags, dim, tags[k]);
- // looping over element types, enlever les déclarations des boucles pr aller plus vite (cfr Boman)
+ // looping over element types
for (std::size_t i = 0; i < elementTypes.size(); ++i)
{
// getting element properties
@@ -1446,9 +1559,7 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
std::vector<double> jacobians, determinants, coords;
gmsh::model::mesh::getJacobians(elementTypes[i],
localCoords, jacobians,
- determinants, coords, tags[k]); //big modif
-
- //std::cout << "entity " << k << ", element type " << i << ",jacobians size: " << jacobians.size() << ", weights size: " << weights.size() << "\n";
+ determinants, coords, tags[k]);
// values and derivatives of the shape functions at the Gauss points of the reference element
// we are here working with derivatives in the reference space, same for all elements
@@ -1462,12 +1573,8 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
basisFunctionsGradient,
numOrientations);
- //std::cout << "entity " << k << ", element tags " << elementTags[i].size() << ",jacobians size: " << jacobians.size() << ",determinants size: " << determinants.size() << ", weights size: " << weights.size() << ", basis functions size: " << basisFunctions.size() << ", basis functions gradient size: " << basisFunctionsGradient.size() << "\n";
-
/*---------computing the K matrix (and f vector, only volumic part for now)------*/
// looping over the elements (computing elemental K matrices)
- //double test1 = omp_get_wtime();
- //Eigen::setNbThreads(1);
#ifdef _MSC_VER
// https://github.com/tesseract-ocr/tesseract/issues/3285
// OpenMP 2.0 standard is enforced - we may use /openmp:llvm and keep size_t
@@ -1475,7 +1582,8 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
for (int el = 0; el < elementTags[i].size(); el++){
#else
#pragma omp parallel for
- for (std::size_t el = 0; el < elementTags[i].size(); el++){
+ for (std::size_t el = 0; el < elementTags[i].size(); el++)
+ {
#endif
// elemental K matrix and f vector initialization
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> K_matrix(2*numNodes,2*numNodes);
@@ -1483,19 +1591,17 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
K_matrix.setZero();
for (int j = 0; j < 2*numNodes; ++j)f_vector[j] = 0.;
- // looping over the Gauss points (formula of slide 17)
- for (std::size_t j = 0; j < weights.size(); ++j){
+ // looping over the Gauss points
+ for (std::size_t gaussIndex = 0; gaussIndex < weights.size(); ++gaussIndex){
// converting jacobian to Eigen format
Eigen::Matrix<double, 2, 2> jacob2D;
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
- jacob(&jacobians[9*j+9*weights.size()*el], 3, 3);
+ jacob(&jacobians[9*gaussIndex+9*weights.size()*el], 3, 3);
// 3D to 2D (plane stress case)
jacob2D(0,0) = jacob(0,0);
jacob2D(1,0) = jacob(1,0);
jacob2D(0,1) = jacob(0,1);
jacob2D(1,1) = jacob(1,1);
-
- //std::cout << "element tag: " << elementTags[i][el] << "jacob2d: " << jacob2D << "\n";
// computing the inverse
Eigen::Matrix<double, 2, 2> jacobinv = jacob2D.inverse();
@@ -1503,41 +1609,107 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
// computing the transpose of the inverse
Eigen::Matrix<double, 2, 2> jacobinvtrans = jacobinv.transpose();
- // computing the B matrix (formula slide 19) and the f vector (also using gauss integration)
+ // computing the B matrix and the f vector (also using gauss integration)
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
-
- computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
-
- //std::cout << "element tag: " << elementTags[i][el] << "B matrix: " << B_matrix << "\n";
+ computeBmatrix(gaussIndex, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
for (int l = 0; l < numNodes; ++l){ // looping over the number of shape functions
- f_vector[2*l] = f_vector[2*l] + rho*basisFunctions[l+numNodes*j]*bx*determinants[j+weights.size()*el]*weights[j]; //pas oublier le dét !
- f_vector[2*l+1] = f_vector[2*l+1] + rho*basisFunctions[l+numNodes*j]*by*determinants[j+weights.size()*el]*weights[j]; // pas oublier le weights[j] non plus !
+ f_vector[2*l] = f_vector[2*l] + rho*basisFunctions[l+numNodes*gaussIndex]*bx*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
+ f_vector[2*l+1] = f_vector[2*l+1] + rho*basisFunctions[l+numNodes*gaussIndex]*by*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
}
- K_matrix = K_matrix + B_matrix.transpose()*H_matrix*B_matrix*determinants[j+weights.size()*el]*weights[j];
- //std::cout << "determinant: " << determinants[j+weights.size()*el] << ", weight: " << weights[j] << "\n";
+ K_matrix = K_matrix + B_matrix.transpose()*H_matrix*B_matrix*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
}
/*------------- ASSEMBLY PROCESS ------------*/
// USE LIST OF TRIPLETS TO PUSH ONLY ONCE -> more rapid and easier to parallelize
- assembleKtriplet(i, el, nodeTags, numNodes, nodeIndexMap, my_k_list[omp_get_thread_num()], K_matrix, f_vector);
- assembleFthread(i, el, nodeTags, numNodes, nodeIndexMap, thread_f_vector[omp_get_thread_num()], f_vector);
+ assembleKtriplet(el, nodeTags[i], numNodes, nodeIndexMap, my_k_list[omp_get_thread_num()], K_matrix);
+ assembleFthread(el, nodeTags[i], numNodes, nodeIndexMap, thread_f_vector[omp_get_thread_num()], f_vector);
+ }
+ }
+ }
+ // gathering all triplet vectors in one single vector
+ std::vector<Eigen::Triplet <double>> k_tripletList = my_k_list[0];
+ if(max_threads > 1)
+ {
+ for (int i = 1; i < max_threads; i++)
+ {
+ k_tripletList.insert(k_tripletList.end(), my_k_list[i].begin(), my_k_list[i].end());
+ }
+ }
+ // assembling the full sparse K matrix
+ full_K_matrix.setFromTriplets(k_tripletList.begin(), k_tripletList.end());
+
+ // re-assembling the full f vector based on the separated thread vectors
+ for(int i = 0; i < max_threads; i++){
+ std::vector<std::pair<int,double>> &tmp_vector = thread_f_vector[i];
+ for(auto &p : tmp_vector) {
+ full_f_vector[p.first] += p.second;
+ }
+ }
+}
+
+// assembles the elemental "K_matrix" into a vector of triplets "k_tripletList" using
+// the node tags provided by "nodeTags" vector (index of the element is "elIndex").
+// The element has "numNodes" nodes. "nodeIndexMap" provides the global nodal index corresponding to a given node tag.
+void assembleKtriplet(const int &elIndex,
+ const std::vector<std::size_t> & nodeTags,
+ const int numNodes, std::map<int, int> & nodeIndexMap,
+ std::vector<Eigen::Triplet <double>> & k_tripletList,
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix)
+{
+ //we go through the lower half of the elemental K_matrix before assembling in the triplet vector
+ for (int j = 0; j < 2*numNodes; ++j){
+ int first_node_index = (int) j/2;
+ int first_node = nodeTags[numNodes*elIndex+first_node_index];
+ int my_row = nodeIndexMap[first_node] + j % 2;
+ for(int l = 0; l <= j; ++l){ // to exploit symmetry
+ int second_node_index = (int) l/2;
+ int second_node = nodeTags[numNodes*elIndex+second_node_index];
+ int my_col = nodeIndexMap[second_node] + l % 2;
- //std::cout << "element tag: " << elementTags[i][el] << "local K matrix: " << K_matrix << "\n";
+ k_tripletList.push_back(Eigen::Triplet<double>(my_row,my_col,K_matrix(j,l)));
+ if(j != l){
+ k_tripletList.push_back(Eigen::Triplet<double>(my_col,my_row,K_matrix(j,l)));
}
- //double test2 = omp_get_wtime() - test1;
- //std::cout << "loop time: " << test2 << "\n";
}
}
}
+// assembles the elemental "f_vector" in a vector of pairs called "thread_f_vector" particular to the parallel thread.
+// the element index in "nodeTags" is "elIndex". The element has "numNodes" nodes and their global nodal indices are given
+// by the "nodeIndexMap".
+void assembleFthread(const int &elIndex,
+ const std::vector<std::size_t> & nodeTags,
+ const int &numNodes, std::map<int, int> & nodeIndexMap,
+ std::vector<std::pair<int,double>> & thread_f_vector,
+ std::vector<double> & f_vector)
+{
+ // we simply go through the elemental f vector
+ int first_node_index;
+ int first_node;
+ int my_row;
+ for (int j = 0; j < 2*numNodes; ++j)
+ {
+ first_node_index = (int) j/2;
+ first_node = nodeTags[numNodes*elIndex+first_node_index];
+ my_row = nodeIndexMap[first_node] + j % 2;
+ thread_f_vector.push_back(std::pair<int, double>(my_row, f_vector[j]));
+ }
+}
+
+// the "full_f_vector" and "full_K_matrix" are reduced in "reduced_f_vector" and "reduced_K_matrix" respectively,
+// taking into account the dirichlet boundary conditions contained in "dir_BC".
+// "new_DOFindices" gathers the new indices of the DOFs when taking the boundary conditions into account.
void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::SparseMatrix<double> & full_K_matrix,
const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC,
Eigen::SparseMatrix<double> & reduced_K_matrix,
- std::vector<double> & reduced_f_vector,
- std::list<Eigen::Triplet <double>> & reduced_k_tripletList)
+ std::vector<double> & reduced_f_vector)
{
+ //the correction term on the f vector comes from the "Finite Element Method" theoretical course (JP Ponthot, ULiege), Chap2 Slide 38/50
+ std::list<Eigen::Triplet <double>> reduced_k_tripletList;
+
+ // initializing the reduced f vector
for (std::size_t j = 0; j < reduced_f_vector.size(); ++j){
reduced_f_vector[j] = 0.;
}
@@ -1547,7 +1719,7 @@ void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::Spars
int new_row;
int new_col;
double k_value;
- for (int col=0; col<full_K_matrix.outerSize(); ++col){
+ for (int col=0; col < full_K_matrix.outerSize(); ++col){
new_col = new_DOFindices[col];
for (Eigen::SparseMatrix<double>::InnerIterator it(full_K_matrix,col); it; ++it){
k_value = it.value();
@@ -1573,12 +1745,42 @@ void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::Spars
reduced_K_matrix.setFromTriplets(reduced_k_tripletList.begin(), reduced_k_tripletList.end());
}
-void computeStressStrainLinearWay(int & dim, std::vector<int> & tags,
- const std::vector<double> &full_d_vector, std::vector<std::size_t> & elems2D,
+
+// computing the elemental-nodal strains and stresses in the FEM domain (dimension "dim" and entity tags "tags"),
+// based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", and the full nodal
+// displacement vector "full_d_vector" and the global nodal index map "nodeIndexMap".
+// "nbelems" is the number of 2D elements contained in the FEM domain.
+// returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag".
+void computeStressStrainLinearWay(int & dim, std::vector<int> & tags, const std::vector<double> &full_d_vector,
const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E,
- std::map<int, int> & nodeIndexMap, std::vector<std::vector<std::vector<double>>> & dataEps,
- std::vector<std::vector<std::vector<double>>> & dataSig)
+ std::map<int, int> & nodeIndexMap, int & nbelems, double & max_VM_stress, int & max_VM_nodeTag)
{
+ std::vector<std::string> names;
+ gmsh::model::list(names);
+
+ int viewTagEpsX = gmsh::view::add("nodal_eps_x");
+ int viewTagEpsY = gmsh::view::add("nodal_eps_y");
+ int viewTagEpsXY = gmsh::view::add("nodal_2eps_xy");
+ int viewTagEpsZ = gmsh::view::add("nodal_eps_z");
+ int viewTagSigX = gmsh::view::add("nodal_sig_x");
+ int viewTagSigY = gmsh::view::add("nodal_sig_y");
+ int viewTagSigXY = gmsh::view::add("nodal_sig_xy");
+ int viewTagSigVM = gmsh::view::add("nodal_sig_VM");
+
+ std::vector<std::vector<std::vector<double>>> dataEps(4);
+ dataEps[0].resize(nbelems); // eps X
+ dataEps[1].resize(nbelems); // eps Y
+ dataEps[2].resize(nbelems); // eps XY
+ dataEps[3].resize(nbelems); // eps Z
+
+ std::vector<std::vector<std::vector<double>>> dataSig(4);
+ dataSig[0].resize(nbelems); // Sig X
+ dataSig[1].resize(nbelems); // Sig Y
+ dataSig[2].resize(nbelems); // Sig XY
+ dataSig[3].resize(nbelems); // Sig VM
+ std::vector<std::size_t> elems2D(nbelems);
+
+
int data_index = 0;
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
@@ -1673,7 +1875,7 @@ void computeStressStrainLinearWay(int & dim, std::vector<int> & tags,
// computing the transpose of the inverse
jacobinvtrans = jacobinv.transpose();
- // computing the B matrix (formula slide 19)
+ // computing the B matrix
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
@@ -1684,18 +1886,46 @@ void computeStressStrainLinearWay(int & dim, std::vector<int> & tags,
dataEps[0][data_index][j] = epsilon(0);
dataEps[1][data_index][j] = epsilon(1);
dataEps[2][data_index][j] = epsilon(2);
- dataEps[3][data_index][j] = -nu/E*(sigma(0)+sigma(1));// formula slide 11/20 elasticity
+ dataEps[3][data_index][j] = -nu/E*(sigma(0)+sigma(1));
dataSig[0][data_index][j] = sigma(0);
dataSig[1][data_index][j] = sigma(1);
dataSig[2][data_index][j] = sigma(2);
dataSig[3][data_index][j] = sqrt(sigma(0)*sigma(0) + sigma(1)*sigma(1) - sigma(0)*sigma(1) + 3*sigma(2)*sigma(2));
+
+ //retrieving maximal nodal von mises stress
+ if(dataSig[3][data_index][j] > max_VM_stress)
+ {
+ max_VM_stress = dataSig[3][data_index][j];
+ max_VM_nodeTag = (int) elnodeTags[ityp][numNodes*i+j];
+ }
}
data_index++;
}
}
}
+
+ gmsh::view::addModelData(viewTagEpsX, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[0], 1);
+ gmsh::view::addModelData(viewTagEpsY, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[1], 1);
+ gmsh::view::addModelData(viewTagEpsXY, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[2], 1);
+ gmsh::view::addModelData(viewTagEpsZ, 0, names[0], "ElementNodeData",
+ elems2D, dataEps[3], 1);
+ gmsh::view::addModelData(viewTagSigX, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[0], 1);
+ gmsh::view::addModelData(viewTagSigY, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[1], 1);
+ gmsh::view::addModelData(viewTagSigXY, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[2], 1);
+ gmsh::view::addModelData(viewTagSigVM, 0, names[0], "ElementNodeData",
+ elems2D, dataSig[3], 1);
}
+// computing the total (internal) strain energy "strain_energy" in the whole FEM domain (dimension "dim" and
+// entity tags "tags"), using the nodal displacements in "full_d_vector". The nodes are indexed using the "nodeIndexMap".
+// The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains.
+// Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule".
void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<int> & tags,
const std::vector<double> &full_d_vector, const Eigen::Matrix<double, 3, 3> &H_matrix,
std::map<int, int> & nodeIndexMap, const std::string & integration_rule)
@@ -1754,7 +1984,6 @@ void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<
Eigen::VectorXd nodal_displacement(2*numNodes); //use eigen vector to use Eigen matric product later
-
for (std::size_t el = 0; el < elementTags[i].size(); el++){
for (int j = 0; j < numNodes; ++j)
@@ -1766,12 +1995,12 @@ void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<
double elementalStrainEnergy = 0;
- // looping over the Gauss points (formula of slide 17)
- for (std::size_t j = 0; j < weights.size(); ++j){
+ // looping over the Gauss points
+ for (std::size_t gaussIndex = 0; gaussIndex < weights.size(); ++gaussIndex){
// converting jacobian to Eigen format
Eigen::Matrix<double, 2, 2> jacob2D;
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
- jacob(&jacobians[9*j+9*weights.size()*el], 3, 3);
+ jacob(&jacobians[9*gaussIndex+9*weights.size()*el], 3, 3);
// 3D to 2D (plane stress case)
jacob2D(0,0) = jacob(0,0);
jacob2D(1,0) = jacob(1,0);
@@ -1784,11 +2013,10 @@ void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<
// computing the transpose of the inverse
Eigen::Matrix<double, 2, 2> jacobinvtrans = jacobinv.transpose();
- // computing the B matrix (formula slide 19) and the f vector (also using gauss integration)
+ // computing the B matrix
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
-
- computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
+ computeBmatrix(gaussIndex, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
epsilon = B_matrix*nodal_displacement;
sigma = H_matrix*epsilon;
@@ -1799,7 +2027,7 @@ void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<
dot_product += epsilon(z)*sigma(z);
}
- elementalStrainEnergy += 0.5*dot_product*determinants[j+weights.size()*el]*weights[j];
+ elementalStrainEnergy += 0.5*dot_product*determinants[gaussIndex+weights.size()*el]*weights[gaussIndex];
}
/*------------- ASSEMBLY PROCESS ------------*/
strainEnergy += elementalStrainEnergy;
diff --git a/srcs/FEM/functionsFEM.hpp b/srcs/FEM/functionsFEM.hpp
index 4572228..1e9df82 100644
--- a/srcs/FEM/functionsFEM.hpp
+++ b/srcs/FEM/functionsFEM.hpp
@@ -1,10 +1,11 @@
#include <iostream>
-#include <vector> // Why should I include this ?
+#include <vector>
#include <map>
#include <Eigen/Dense>
#include <list>
-#include <Eigen/Sparse> // Eigen library for sparse matrices
+#include <Eigen/Sparse>
+// data structure used for the non-linear solver "solverFEMnonLinear"
struct elementData{
// will store all the kinematics data
std::vector<size_t> nodes; // will contain the tags of the nodes
@@ -34,47 +35,253 @@ struct elementData{
Eigen::Matrix<double, Eigen::Dynamic, 1> fl; //local current (iterative) forces under vector form
- Eigen::Matrix<double, Eigen::Dynamic, 1> A;
- Eigen::Matrix<double, 1, Eigen::Dynamic> G;
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> P;
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> E;
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B;
+ Eigen::Matrix<double, Eigen::Dynamic, 1> A; // see report for exact definition
+ Eigen::Matrix<double, 1, Eigen::Dynamic> G; // see report for exact definition
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> P; // see report for exact definition
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> E; // see report for exact definition
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> C; // see report for exact definition
};
+/*------------ "SOLVER" functions -----------------*/
+
+// solves the non-linear (large displacements, no large deformations) elastic problem.
+// if coupled to the BEM solver, "electrostaticPressure" is the pressure computed from the electric field
+// on the BEM_FEM boundary (mapping element tag -> pressure). "boundaryDisplacementMap" is filled for all
+// nodes on the BEM_FEM_boundary as (node tag -> (u,v) displacement)
+// "nbViews" is the current number of views already allocated in the gmsh window.
+// "postProcessing" boolean activates the post processing operations (must only be computed once at the end of the
+// coupled solver). "iteration" denotes the current iteration of the two-way coupled solver (starting at 1).
+// "viewTagU" and "viewTagF" correspond to the view tags of the nodal displacements and the nodal forces in the gmsh window.
+// when "untangleMesh" is activated, the nodal coordinates of the nodes in the FEM_domain are updated according to
+// their final nodal displacements.
+void solverFEMnonLinear(std::map<int, double> &electrostaticPressure,
+ std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int &nbViews, bool postProcessing,
+ const int iteration, const int viewTagU, const int viewTagF, bool untangleMesh);
+
+// solves the linear elastic problem in the FEM domain.
+// "electrostaticPressure" contains the pressure resulting from the electric field computation on the
+// BEM_FEM_boundary. "nbViews" is the current number of views already allocated in the gmsh window.
+void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews);
+
+/*------------ More specific functions called inside the "SOLVER" functions -----------------*/
+
+// returns the boolean non_linear_solver parameter based on what is set in the .geo file,
+// add SetNumber("Non_linear_solver", 0); in .geo file to use linear solver.
bool getNonLinearParameter();
-Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu);
-void assembleFlocal(const int &i, const int &el, const std::vector<std::vector<std::size_t>> & nodeTags, const int &numNodes, std::map<int, int> & nodeIndexMap, std::vector<double> & full_f_vector, std::vector<double> & f_vector);
-void computeBmatrix(const int &j, const int &numNodes, const std::vector<double> & basisFunctionsGradient, const std::vector<double> & determinants, const Eigen::Matrix<double, 2, 2> & jacobinvtrans, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & B_matrix);
-void getPhysicalProperties(const std::vector<std::string> & keys, double & E, double & nu, double & rho, double & bx, double & by);
-void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, const std::string & integration_rule, int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap, std::vector<double> & full_f_vector, std::map<int, elementData> & FEM_kinematicsMap);
-void computeNeumannBC(const int &numNodes, const int &elTyp, const std::vector<std::vector<std::size_t>> & elementTags, const double &tx, const double &ty, std::map<int, int> & nodeIndexMap, const std::vector<double> & weights, const std::vector<double> & basisFunctions, const std::vector<double> & determinants, const std::vector<std::vector<std::size_t>> & nodeTags, std::vector<double> & full_f_vector);
-void includeNeumannBC(const std::vector<std::string> & keys, const std::string & integration_rule, std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups, std::vector<double> & full_f_vector);
-void fillDOFindicesDirBC(const std::vector<std::string> & keys, std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups, std::vector<int> & new_DOFindices, std::vector<double> & dir_BC);
-void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups, std::map<int, int> & nodeIndexMap, std::vector<double> & nodal_forces, const std::vector<std::string> & keys);
-void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector, std::vector<std::string> & names);
-void computeStressStrainElNodal(int & dim, std::vector<int> & tags, std::vector<std::size_t> & elems2D, const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E, std::map<int, elementData> & FEM_kinematicsMap, std::vector<std::vector<std::vector<double>>> & dataEps, std::vector<std::vector<std::vector<double>>> & dataSig, double & max_VM_stress, int & max_VM_nodeTag);
-void applyElecPressure(const std::string & integration_rule, std::map<int,std::pair<double,double>> & nodeCoordsMap, std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups, std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector);
-int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, elementData> & FEM_kinematicsMap, std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags);
-void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector, std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags);
+
+// fills the "FEM_nodeTags" vector with the nodes in the FEM domain, fills the "FEM_nodeCoordsMap" with the coordinates (x,y)
+// corresponding to each nodeTag, looking at entity tags "tags" of dimension "dim" corresponding to the FEM domain.
+void RetrieveFEMNodeTagsAndCoords(std::vector<std::size_t> & FEM_nodeTags,
+ std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap,
+ int & dim, std::vector<int> & tags);
+
+// fills the "FEM_elementTags" vector with all element tags of the elements of dimension "dim" (in practice: 2)
+// in the entities associated to the entity tags "tags" (corresponding to the FEM domain).
+// associates an initialized elementData structure to each element tag through the "FEM_kinematicsMap" map,
+// using the "FEM_nodeCoordsMap" map gathering the coordinates of each node in the FEM domain.
+int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags);
+
+// updates the kinematics matrices of one element "element" based on its updated nodal positions and rotation angle.
+// if "computeFl" is set as true, also computes the local nodal forces vector "fl = Kl*pl", only OK if
+// the local stiffness matrix Kl has already been initialized, using "fillElementalKandF".
void updateKinMatrices(elementData* element, bool computeFl);
-void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elementData> & FEM_kinematicsMap, std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, std::vector<std::size_t> & FEM_elementTags);
-double computeResidualNorm(std::vector<double> & global_f_vector, std::vector<double> & f_app_vector, Eigen::Matrix<double, Eigen::Dynamic, 1> & residual);
-void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::map<int, elementData> & FEM_kinematicsMap, std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, std::vector<std::size_t> & FEM_elementTags);
-void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<double> & global_p_vector, std::vector<int> & new_DOFindices);
-void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map<int, elementData> & FEM_kinematicsMap, std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags);
-void retrieveBoundaryDisplacements(std::map<int,std::pair<double,double>> & boundaryDisplacementMap, const std::vector<double> & full_d_vector, std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups);
-void displayTime(double &start, double &end, std::string &function_name, bool &display_time);
-void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & full_d_vector, std::vector<double> & final_nodal_forces);
-void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & tags, std::map<int, elementData> & FEM_kinematicsMap, const Eigen::Matrix<double, 3, 3> &H_matrix, const std::string & integration_rule);
-void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int &nbViews, bool postProcessing, const int iteration, const int viewTagU, const int viewTagF, bool untangleMesh);
+// gets physical properties of FEM_domain: Young's modulus "E", Poisson's ratio "nu",
+// density "rho", volumic force along x "bx", volumic force along y "by".
+void getPhysicalProperties(double & E, double & nu, double & rho, double & bx, double & by);
+
+// computes the Hooke's matrix in plane stress "H_matrix" based on the Young's modulus "E" and Poisson's ratio "nu".
+Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu);
+
+// looping over every element in the FEM_domain (dimension "dim" and entity tags "tags"), filling "full_f_vector",
+// with the local volumic forces computed based on the volumic force ("bx","by") and density "rho" with the help of
+// "nodeIndexMap" mapping each nodeTag to its global index.
+// filling the local stiffness matrices stored in the elementData structure contained in "FEM_kinematicsMap".
+// numerical integration performed using Gauss integration based on "integration_rule".
+void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by,
+ const std::string & integration_rule,
+ int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap,
+ std::vector<double> & full_f_vector,
+ std::map<int, elementData> & FEM_kinematicsMap);
+
+// computing the strain displacement "B_matrix" at the "gaussIndex"-th gaussian point of the element
+// consisting of "numNodes" nodes, based on the "basisFunctionsGradient", "determinants" and the inverse of the jacobian
+// matrix "jacobinvtrans".
+void computeBmatrix(const int &gaussIndex, const int &numNodes, const std::vector<double> & basisFunctionsGradient,
+ const std::vector<double> & determinants,
+ const Eigen::Matrix<double, 2, 2> & jacobinvtrans,
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & B_matrix);
+
+// assembles the local "f_vector" into the full f vector, "elIndex" is the element index with respect to "nodeTags"
+// having "numNodes" nodes, "nodeTags" gathering all the nodeTags of the current element type in the current entity.
+// "nodeIndexMap" associates a global index to each nodetag.
+void assembleFlocal(const int &elIndex,
+ const std::vector<std::size_t> & nodeTags,
+ const int &numNodes, std::map<int, int> & nodeIndexMap,
+ std::vector<double> & full_f_vector,
+ std::vector<double> & f_vector);
+
+// fills the "full_f_vector" with the surfacic contribution of the neumann boundary conditions based on the
+// "nodeIndexMap" associating a global index to each node tag, the local integration of the surfacic forces is
+// performed using Gauss integration based on "integration rule".
+// "BCkeys" gathers the keys associated to boundary conditions in the .geo file (using SetNumber(...)).
+// "groups" associates a dimension and a tag to each physical group.
+void includeNeumannBC(const std::vector<std::string> & BCkeys, const std::string & integration_rule,
+ std::map<int, int> & nodeIndexMap,
+ std::map<std::string, std::pair<int, int>> & groups,
+ std::vector<double> & full_f_vector);
+
+// for the coupled solver, fills the "full_f_vector" with the electrostatic pressure contained in the
+// "electrostaticPressure" map for different element tags, using Gauss integration "integration_rule", the "nodeIndexMap"
+// associating a global nodal index to each nodeTag, "nodeCoordsMap" gathering the coordinates of each node and "groups"
+// containing the different physical groups of the domain.
+void applyElecPressure(const std::string & integration_rule, std::map<int,std::pair<double,double>> & nodeCoordsMap,
+ std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups,
+ std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector);
+
+// looping over the boundary conditions imposed in the .geo file and gathered in "BCkeys" related to the
+// physical groups "groups", the goal is to fill the two vectors "new_DOFindices" and "dir_BC":
+// "new_DOFindices" vector will contain the new indices of the K matrix after the "number_removed_lines" rows/columns
+// corresponding to dirichlet boundary conditions have been removed. A removed index will correspond to "new_DOFindices"= -1.
+// dir_BC[index] will contain the dirichlet boundary condition prescribed to the node whose global index is retrieved using
+// "nodeIndexMap".
+void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys,
+ std::map<int, int> & nodeIndexMap,
+ std::map<std::string, std::pair<int, int>> & groups,
+ std::vector<int> & new_DOFindices,
+ std::vector<double> & dir_BC, int & number_removed_lines);
+
+// updating the kinematics of the "FEM_elementTags" contained in the "FEM_kinematicsMap",
+// based on the "full_d_vector" gathering the nodal displacements associated to the global
+// node indices contained in "nodeIndexMap".
+void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector,
+ std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags);
+
+// retrieving the global nodal forces based on local nodal forces of the
+// "FEM_elementTags" contained in the "FEM_kinematicsMap",
+// each node Tag corresponds to a global index given by "nodeIndexMap".
+void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags);
+
+// assembles the "global_f_vector" based on the local nodal forces vector of the
+// "FEM_elementTags" contained in the "FEM_kinematicsMap",
+// without taking into account the DOFs fixed by a dirichlet boundary condition,
+// the nodes are renumbered using "new_DOFindices".
+// each node Tag corresponds to a global index given by "nodeIndexMap".
+void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices,
+ std::vector<std::size_t> & FEM_elementTags);
+
+// assembling the global tangent stiffness matrix based on the local elemental stiffness matrices contained
+// in the "FEM_kinematicsMap" associated to each "FEM_elementTags".
+// "new_DOFindices" maps the global initial node index (from "nodeIndexMap") to the index
+// obtained without considering the fixed DOFs.
+void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::map<int, elementData> & FEM_kinematicsMap,
+ std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices,
+ std::vector<std::size_t> & FEM_elementTags);
+
+// updating the "full_d_vector" based on the "global_p_vector", using the global and reduced index mapping:
+// "new_DOFindices".
+void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<double> & global_p_vector,
+ std::vector<int> & new_DOFindices);
+
+// retrieving the nodal displacements on the BEM_FEM_boundary in order to communicate them to the
+// BEM solver in the two-way coupled solver. Using the "groups" to find the physical group associated
+// to the boundary, all nodes of the boundary are associated through the "boundaryDisplacementMap" their
+// nodal displacement (u,v) retrieved from "full_d_vector" using the indexation of "nodeIndexMap".
+void retrieveBoundaryDisplacements(std::map<int,std::pair<double,double>> & boundaryDisplacementMap,
+ const std::vector<double> & full_d_vector, std::map<int, int> & nodeIndexMap,
+ std::map<std::string, std::pair<int, int>> & groups);
+
+// computing the reaction forces on all the physical groups "groups" for which a dirichlet boundary
+// condition has been prescribed (using the "BCkeys" to retrieve the groups from the .geo file),
+// based on the "final_nodal_forces" and the "nodeIndexMap" to retrieve the global index associated
+// to one given node tag.
+void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups,
+ std::map<int, int> & nodeIndexMap, std::vector<double> & final_nodal_forces,
+ const std::vector<std::string> & BCkeys);
+
+// plotting the nodal displacement values associated to "FEM_nodeTags" in the "names[0]" model,
+// based on the "full_d_vector".
+void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector,
+ std::vector<std::string> & names);
+
+// computing the elemental-nodal strains and stresses in the FEM domain (dimension "dim" and entity tags "tags"),
+// based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", looping over the
+// "nbelems" elements contained in "FEM_kinematicsMap".
+// returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag".
+void computeStressStrainElNodal(int & dim, std::vector<int> & tags, const Eigen::Matrix<double, 3, 3> &H_matrix,
+ const double nu, const double E, std::map<int, elementData> & FEM_kinematicsMap,
+ int & nbelems, double & max_VM_stress, int & max_VM_nodeTag);
+
+// computing the work done by external forces, as the nodal forces are energy consistent, it
+// is sufficient to compute the "externalWork" as the dot product of the nodal displacements
+// "full_d_vector" and the nodal forces "final_nodal_forces".
+void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & full_d_vector,
+ std::vector<double> & final_nodal_forces);
+
+// computing the total (internal) strain energy "strain_energy" in the whole FEM domain (dimension "dim" and
+// entity tags "tags"), using the elementData stored in "FEM_kinematicsMap" gathering the informations for all
+// FEM elements. The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains.
+// Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule".
+void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & tags,
+ std::map<int, elementData> & FEM_kinematicsMap, const Eigen::Matrix<double, 3, 3> &H_matrix,
+ const std::string & integration_rule);
+
+//displays the time of the desired code section "function_name", if "display_time" is set as true.
+void displayTime(double &start, double &end, std::string &function_name, bool &display_time);
/*-------------- FUNCTIONS SPECIFIC TO LINEAR FEM ------------------*/
-void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews);
-void assembleKtriplet(const int i, const int el, const std::vector<std::vector<std::size_t>> & nodeTags, const int numNodes, std::map<int, int> & nodeIndexMap, std::vector<Eigen::Triplet <double>> & k_tripletList, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix, std::vector<double> & f_vector);
-void assembleFthread(const int &i, const int &el, const std::vector<std::vector<std::size_t>> & nodeTags, const int &numNodes, std::map<int, int> & nodeIndexMap, std::vector<std::pair<int,double>> & thread_f_vector, std::vector<double> & f_vector);
-void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, const std::string & integration_rule, int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap, std::vector<std::vector<Eigen::Triplet <double>>> & my_k_list, std::vector<std::vector<std::pair<int,double>>> & thread_f_vector);
-void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::SparseMatrix<double> & full_K_matrix, const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC, Eigen::SparseMatrix<double> & reduced_K_matrix, std::vector<double> & reduced_f_vector, std::list<Eigen::Triplet <double>> & reduced_k_tripletList);
-void computeStressStrainLinearWay(int & dim, std::vector<int> & tags, const std::vector<double> &full_d_vector, std::vector<std::size_t> & elems2D, const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E, std::map<int, int> & nodeIndexMap, std::vector<std::vector<std::vector<double>>> & dataEps, std::vector<std::vector<std::vector<double>>> & dataSig);
-void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<int> & tags, const std::vector<double> &full_d_vector, const Eigen::Matrix<double, 3, 3> &H_matrix, std::map<int, int> & nodeIndexMap, const std::string & integration_rule);
\ No newline at end of file
+// assembling the full stiffness matrix "full_K_matrix" using first a triplet vector, assembling the volumic part
+// of the nodal forces "full_f_vector", in the FEM_domain (dimension "dim" and entity tag "tags") using the global
+// node indexation defined by "nodeIndexMap". "H_matrix" is the plane stress Hooke's matrix, "rho" the density of the
+// material, "bx" the volumic force along x, "by" the volumic force along y.
+// Numerical integration using gauss integration "integration_rule".
+void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by,
+ const std::string & integration_rule, int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap,
+ Eigen::SparseMatrix<double> & full_K_matrix, std::vector<double> & full_f_vector);
+
+// assembles the elemental "K_matrix" into a vector of triplets "k_tripletList" using
+// the node tags provided by "nodeTags" vector (index of the element is "elIndex").
+// The element has "numNodes" nodes. "nodeIndexMap" provides the global nodal index corresponding to a given node tag.
+void assembleKtriplet(const int &elIndex,
+ const std::vector<std::size_t> & nodeTags,
+ const int numNodes, std::map<int, int> & nodeIndexMap,
+ std::vector<Eigen::Triplet <double>> & k_tripletList,
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix);
+
+// assembles the elemental "f_vector" in a vector of pairs called "thread_f_vector" particular to the parallel thread.
+// the element index in "nodeTags" is "elIndex". The element has "numNodes" nodes and their global nodal indices are given
+// by the "nodeIndexMap".
+void assembleFthread(const int &elIndex,
+ const std::vector<std::size_t> & nodeTags,
+ const int &numNodes, std::map<int, int> & nodeIndexMap,
+ std::vector<std::pair<int,double>> & thread_f_vector,
+ std::vector<double> & f_vector);
+
+// the "full_f_vector" and "full_K_matrix" are reduced in "reduced_f_vector" and "reduced_K_matrix" respectively,
+// taking into account the dirichlet boundary conditions contained in "dir_BC".
+// "new_DOFindices" gathers the new indices of the DOFs when taking the boundary conditions into account.
+void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::SparseMatrix<double> & full_K_matrix,
+ const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC,
+ Eigen::SparseMatrix<double> & reduced_K_matrix,
+ std::vector<double> & reduced_f_vector);
+
+// computing the elemental-nodal strains and stresses in the FEM domain (dimension "dim" and entity tags "tags"),
+// based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", and the full nodal
+// displacement vector "full_d_vector" and the global nodal index map "nodeIndexMap".
+// "nbelems" is the number of 2D elements contained in the FEM domain.
+// returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag".
+void computeStressStrainLinearWay(int & dim, std::vector<int> & tags, const std::vector<double> &full_d_vector,
+ const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E,
+ std::map<int, int> & nodeIndexMap, int & nbelems, double & max_VM_stress, int & max_VM_nodeTag);
+
+// computing the total (internal) strain energy "strain_energy" in the whole FEM domain (dimension "dim" and
+// entity tags "tags"), using the nodal displacements in "full_d_vector". The nodes are indexed using the "nodeIndexMap".
+// The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains.
+// Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule".
+void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<int> & tags,
+ const std::vector<double> &full_d_vector, const Eigen::Matrix<double, 3, 3> &H_matrix,
+ std::map<int, int> & nodeIndexMap, const std::string & integration_rule);
\ No newline at end of file
diff --git a/srcs/FEM/large_rotation_validation.geo b/srcs/FEM/large_rotation_validation.geo
index f2bcf80..113db4e 100644
--- a/srcs/FEM/large_rotation_validation.geo
+++ b/srcs/FEM/large_rotation_validation.geo
@@ -2,7 +2,7 @@
h = 1;
H = 10;
-n = 48;
+n = 16;
Point(1) = {0, 0, 0, 0.5};
Point(2) = {0.9*H, 0, 0, 0.5};
diff --git a/srcs/FEM/mainFEM.cpp b/srcs/FEM/mainFEM.cpp
index 02b644c..3883e6b 100644
--- a/srcs/FEM/mainFEM.cpp
+++ b/srcs/FEM/mainFEM.cpp
@@ -4,6 +4,8 @@
#include <iostream>
#include <omp.h>
+// this program implements a FEM solver, either linear or non-linear (by taking large displacements into account)
+// to choose the linear solver, please add: SetNumber("Non_linear_solver",0); in the .geo file.
int main(int argc, char **argv)
{
if (argc < 2)
@@ -28,18 +30,19 @@ int main(int argc, char **argv)
// determine if non linear solver must be set
bool nonLinearSolver = getNonLinearParameter();
+ // maps the elementTags to a corresponding electrostaticPressure (only useful for coupled solver)
std::map<int, double> electrostaticPressure;
- int nbViews = 0;
+ int nbViews = 0; // number of views passed to the BEM solver (only useful for coupled solver)
if(nonLinearSolver)
{
- std::map<int,std::pair<double,double>> boundaryDisplacementMap;
+ std::map<int,std::pair<double,double>> boundaryDisplacementMap; //(only useful for coupled solver)
- int viewTagU = gmsh::view::add("u"); // Ã modifier plus tard
- int viewTagF = gmsh::view::add("f");
+ int viewTagU = gmsh::view::add("u"); //viewtag of the displacement view
+ int viewTagF = gmsh::view::add("f"); //viewtag of the nodal forces view
- solverFEMnonLinear(electrostaticPressure, boundaryDisplacementMap, nbViews, true, 1, viewTagU, viewTagF, false); //iteration randomly set to 1
+ solverFEMnonLinear(electrostaticPressure, boundaryDisplacementMap, nbViews, true, 1, viewTagU, viewTagF, false);
}
else
{
diff --git a/srcs/FEM/solverFEM.cpp b/srcs/FEM/solverFEM.cpp
index 05eb54d..b27c1c8 100644
--- a/srcs/FEM/solverFEM.cpp
+++ b/srcs/FEM/solverFEM.cpp
@@ -1,23 +1,35 @@
#include "functionsFEM.hpp"
#ifdef _MSC_VER
-#include <gmsh.h_cwrap> // gmsh main header
+#include <gmsh.h_cwrap>
#else
-#include <gmsh.h> // gmsh main header
+#include <gmsh.h>
#endif
-#include <iostream> // for std::cout
+#include <iostream>
#include <map>
-#include <sstream> // for std::stringstream
-#include <Eigen/Dense> // Eigen library
-#include <Eigen/Sparse> // Eigen library for sparse matrices
-#include <Eigen/SparseCholesky> // solving sparse linear systems
+#include <sstream>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/SparseCholesky>
#include <Eigen/Core>
#include <cmath>
-#include <algorithm> // to sort and merge the FEM node vectors
+#include <algorithm>
#include <omp.h>
-void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int &nbViews, bool postProcessing, const int iteration, const int viewTagU, const int viewTagF, bool untangleMesh)
+// solves the non-linear (large displacements, no large deformations) elastic problem.
+// if coupled to the BEM solver, "electrostaticPressure" is the pressure computed from the electric field
+// on the BEM_FEM boundary (mapping element tag -> pressure). "boundaryDisplacementMap" is filled for all
+// nodes on the BEM_FEM_boundary as (node tag -> (u,v) displacement)
+// "nbViews" is the current number of views already allocated in the gmsh window.
+// "postProcessing" boolean activates the post processing operations (must only be computed once at the end of the
+// coupled solver). "iteration" denotes the current iteration of the two-way coupled solver (starting at 1).
+// "viewTagU" and "viewTagF" correspond to the view tags of the nodal displacements and the nodal forces in the gmsh window.
+// when "untangleMesh" is activated, the nodal coordinates of the nodes in the FEM_domain are updated according to
+// their final nodal displacements.
+void solverFEMnonLinear(std::map<int, double> &electrostaticPressure,
+ std::map<int,std::pair<double,double>> & boundaryDisplacementMap, int &nbViews, bool postProcessing,
+ const int iteration, const int viewTagU, const int viewTagF, bool untangleMesh)
{
/*--------------INITIALIZATION OF ITERATIVE ALGORITHM----------------*/
/* FIRST STEP:
@@ -26,18 +38,15 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
- initialize the full d vector with Dirichlet boundary conditions and set pg = 0
- important remark: in pg we focus only on the free DOFs (we do not focus on the DOFs fixed by dirichlet BCs)*/
- bool display_time = false;
- bool validation = false;
+ bool display_time = false; //displays time of every different step of the algorithm
+ bool validation = false; //retrieves displacement of point A for large_rotation_validation.geo
double start = omp_get_wtime();
- //int max_threads = omp_get_max_threads();
-
/*--------Integration rule used in the code----------*/
- // could be useful to pass it as an argument in the .geo file
std::string integration_rule = "CompositeGauss4"; // tested with other methods, OK too.
- /*--------get physical groups--------------*/
+ /*--------get all physical groups--------------*/
std::map<std::string, std::pair<int, int>> groups;
gmsh::vectorpair dimTags;
gmsh::model::getPhysicalGroups(dimTags);
@@ -51,46 +60,16 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
/*--------get nodes related to FEM-----------------*/
- std::string groupname = "FEM_domain"; // we consider the domain as containing the FEM 2d model
+ std::string groupname = "FEM_domain"; // we consider the FEM_domain as containing the FEM 2d model
int dim = groups[groupname].first;
int tag = groups[groupname].second;
+ // Getting entities of the domain to fill FEM_nodeTags
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
- //loop over the entities of the domain to retrieve all nodeTags by taking care not to take some nodes twice
- //using sort and set_union methods + retrieving the node coordinates and storing them in a map
std::vector<std::size_t> FEM_nodeTags;
- std::vector<double> FEM_nodecoord;
- std::vector<double> FEM_nodeparametricCoord;
std::map<int,std::pair<double,double>> FEM_nodeCoordsMap; // nodeTag -> (x,y)
- // Getting entities of the domain to fill FEM_nodeTags, en faire une fonction ?
- std::vector<int> tags;
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
- std::vector<std::vector<std::size_t>> tmp_nodeTags(tags.size()); // will store the nodeTags associated to one given entity of the domain
- gmsh::model::mesh::getNodes(tmp_nodeTags[0], FEM_nodecoord, FEM_nodeparametricCoord, dim, tags[0], true);
- for(std::size_t j = 0; j < tmp_nodeTags[0].size(); j++)
- {
- FEM_nodeCoordsMap[tmp_nodeTags[0][j]] = std::pair<double, double>(FEM_nodecoord[3*j], FEM_nodecoord[3*j+1]);
- }
- std::sort(tmp_nodeTags[0].begin(), tmp_nodeTags[0].end()); // need to sort the vectors in order to merge them afterwards
- for (std::size_t i = 1; i < tags.size(); i++)
- { // in the case the domain contains multiple entities
- gmsh::model::mesh::getNodes(tmp_nodeTags[i], FEM_nodecoord, FEM_nodeparametricCoord, dim, tags[i], true);
- for(std::size_t j = 0; j < tmp_nodeTags[i].size(); j++)
- {
- FEM_nodeCoordsMap[tmp_nodeTags[i][j]] = std::pair<double, double>(FEM_nodecoord[3*j], FEM_nodecoord[3*j+1]);
- }
- std::sort(tmp_nodeTags[i].begin(), tmp_nodeTags[i].end());
- std::set_union(tmp_nodeTags[i-1].begin(), tmp_nodeTags[i-1].end(),
- tmp_nodeTags[i].begin(), tmp_nodeTags[i].end(), std::back_inserter(FEM_nodeTags)); // works only for 2 domains -> see little trick just below
-
- }
- if(tags.size() == 1){ // if the domain contains one single entity
- FEM_nodeTags = tmp_nodeTags[0];
- }
- else{ // little trick for more than 2 domains
- std::sort(FEM_nodeTags.begin(), FEM_nodeTags.end());
- auto last = std::unique(FEM_nodeTags.begin(), FEM_nodeTags.end());
- FEM_nodeTags.erase(last, FEM_nodeTags.end());
- }
+ RetrieveFEMNodeTagsAndCoords(FEM_nodeTags, FEM_nodeCoordsMap, dim, tags);
/*-----get elements related to FEM--------------*/
std::vector<std::size_t> FEM_elementTags;
@@ -113,15 +92,12 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
/*-----------get physical properties of domain----------------------*/
- std::vector<std::string> keys;
- gmsh::onelab::getNames(keys, "(Volumic Forces|Materials).+");
double E = 0; // Young modulus
double nu = 0; // Poisson ratio
double rho = 0; // volumic mass
double bx = 0; // volumic force along x axis
double by = 0; // volumic force along y axis
-
- getPhysicalProperties(keys, E, nu, rho, bx, by);
+ getPhysicalProperties(E, nu, rho, bx, by);
/*-----------------H matrix involved in stress computation------------------*/
Eigen::Matrix<double, 3, 3> H_matrix = computeHmatrix(E, nu);
@@ -137,7 +113,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
full_f_vector[i] = 0;
}
- // Looping over entities and compute local K matrices and volumic f vector
+ // Looping over entities of FEM_dimain and computing local K matrices and volumic f vector
fillElementalKandF(H_matrix, rho, bx, by, integration_rule, dim, tags, nodeIndexMap, full_f_vector, FEM_kinematicsMap);
double start3 = omp_get_wtime();
@@ -147,12 +123,13 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
/*-----------------BOUNDARY CONDITIONS------------*/
/*-----------get group_names related to BC's----------------------*/
- gmsh::onelab::getNames(keys, "(Boundary Conditions).+");
+ std::vector<std::string> BCkeys;
+ gmsh::onelab::getNames(BCkeys, "(Boundary Conditions).+");
/*-----------first focus on neumann BC's (surface traction) into full_f_vector------------*/
- includeNeumannBC(keys, integration_rule, nodeIndexMap, groups, full_f_vector);
+ includeNeumannBC(BCkeys, integration_rule, nodeIndexMap, groups, full_f_vector);
/* taking ELECTROSTATIC PRESSURE into account */
- if(electrostaticPressure.size() != 0)
+ if(electrostaticPressure.size() != 0) // size == 0 if FEM solver used alone
applyElecPressure(integration_rule, FEM_nodeCoordsMap, nodeIndexMap, groups, electrostaticPressure, full_f_vector);
double start4 = omp_get_wtime();
step_name = "Neumann BC and elec pressure";
@@ -162,26 +139,18 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
// new_DOFindices vector will contain the new indices of the K matrix after the rows/columns corresponding
// to dirichlet boundary conditions have been removed
// IN THIS IMPLEMENTATION WE DECIDE TO REMOVE EVERY ROW/COL ASSOCIATED TO DIR. BC's
- // non-homo. dir BC will induce a correction in the RHS term f of the other rows, while homogeneous not
+ // non-homo. dir BC will induce a correction in the RHS term f of the other rows, while homogeneous not,
// if the given index has been removed, it will be set to -1 in new_DOFindices
std::vector<int> new_DOFindices(2*FEM_nodeTags.size());
// dir_BC[index] will contain the dirichlet boundary condition prescribed to this node
std::vector<double> dir_BC(2*FEM_nodeTags.size());
- // Fill the two vectors defined just above
- fillDOFindicesDirBC(keys, nodeIndexMap, groups, new_DOFindices, dir_BC);
-
+ // Fills the two vectors defined just above
int number_removed_lines = 0;
- for (std::size_t j = 0; j < new_DOFindices.size(); ++j){
- if(new_DOFindices[j] < 0){ // line has to be removed
- number_removed_lines++;
- }
- else{
- new_DOFindices[j] -= number_removed_lines;
- }
- }
+ fillDOFindicesDirBC(BCkeys, nodeIndexMap, groups, new_DOFindices, dir_BC, number_removed_lines);
+ // initializing the full_d_vector with the dirichlet boundary conditions
std::vector<double> full_d_vector(2*FEM_nodeTags.size());
for (std::size_t j = 0; j < full_d_vector.size(); ++j){
if(new_DOFindices[j] < 0){ // line was previously removed due to dirichlet boundary conditions
@@ -191,7 +160,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
full_d_vector[j] = 0;
}
}
- // initializing the applied forces vector (will remain constant) and the global p (p_g) vector
+ // initializing the applied forces vector (will remain constant) and the global p vector (p_g)
std::vector<double> f_app_vector(2*FEM_nodeTags.size() - number_removed_lines);
std::vector<double> global_p_vector(2*FEM_nodeTags.size() - number_removed_lines);
int nb_free_DOFs = 0;
@@ -205,22 +174,11 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
}
- // computing the norm of the applied forces vector
- double f_app_norm = 0;
- for(size_t i = 0; i < f_app_vector.size(); i++)
- {
- f_app_norm += f_app_vector[i]*f_app_vector[i];
- }
- f_app_norm = sqrt(f_app_norm);
- //std::cout << "f app norm: " << f_app_norm << "\n";
-
- std::vector<double> global_f_vector(nb_free_DOFs);
-
+ // useful during the iterative process
+ std::vector<double> global_f_vector(nb_free_DOFs); // f_g
Eigen::Matrix<double, Eigen::Dynamic, 1> residual(nb_free_DOFs, 1);
- double residualTolerance = 1e-12; // empirical
-
- // implementation of second stopping criterion
+ // useful for the stopping criterion
std::vector<double> previousTotalNodalForces(2*FEM_nodeTags.size());
std::vector<double> currentTotalNodalForces(2*FEM_nodeTags.size());
std::vector<double> relativeDifference(2*FEM_nodeTags.size());
@@ -232,6 +190,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
previousTotalNodalForces[i] = 0;
currentTotalNodalForces[i] = 0;
}
+ double residualTolerance = 1e-12; // empirical
double start5 = omp_get_wtime();
step_name = "further initialization of global variables for the iterative algorithm";
@@ -239,7 +198,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
/*---------ITERATIVE PART---------------*/
- int max_number_of_steps = 200; // purely arbitrary
+ int max_number_of_steps = 200; // can be tuned
int step = 0;
double num_diverged_steps = 0; // stores the number of times the residual has increased between two steps, after 5 times,
@@ -289,11 +248,11 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
if(currentMaxNodalForce > 0)
relativeForces = sqrt(relativeForces/(2*FEM_nodeTags.size())/currentMaxNodalForce); // norm made relative
- else
- relativeForces = 1e17;
+ else //to avoid division by zero
+ relativeForces = 1e17;
std::cout << "step: " << step << ", relativeForces: " << relativeForces << "\n";
- if(previousRelativeForces < relativeForces)
+ if(previousRelativeForces < relativeForces) //the current step induces an increase of nodal difference
{
num_diverged_steps++;
if(num_diverged_steps == 5)
@@ -304,7 +263,6 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
}
-
/*--------COMPUTING AND ASSEMBLING THE GLOBAL F VECTOR BASED ON PREVIOUS ITERATION-----------*/
for(int i = 0; i < nb_free_DOFs; i++)
{
@@ -312,8 +270,11 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
assembleGlobalF(global_f_vector, FEM_kinematicsMap, nodeIndexMap, new_DOFindices, FEM_elementTags);
- //double residualNorm = computeResidualNorm(global_f_vector, f_app_vector, residual);
- computeResidualNorm(global_f_vector, f_app_vector, residual);
+ // computing the residual
+ for(size_t i = 0; i < global_f_vector.size(); i++)
+ {
+ residual(i,0) = global_f_vector[i] - f_app_vector[i];
+ }
double end_assembleglobalf = omp_get_wtime();
step_name = "assemble global f and residual norm";
@@ -337,7 +298,8 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
/*-----------STEP 5: SOLVING THE SYSTEM----------*/
Eigen::Matrix<double, Eigen::Dynamic, 1> global_delta_p(nb_free_DOFs, 1);
- // solve Kg*Deltapg=-res
+ // solve Kg*Deltapg=-res, using a direct sparse Cholesky factorization,
+ // efficient for sparse positive semi-definite matrix
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
solver.compute(global_tangent_matrix);
@@ -374,34 +336,13 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
std::vector<double> final_nodal_forces(2*FEM_nodeTags.size());
retrieveTotalNodalForces(final_nodal_forces, FEM_kinematicsMap, nodeIndexMap, FEM_elementTags);
- /*--------------RETRIEVING NODAL DISPLACEMENTS ON BEM-FEM BOUNDARY------------*/
- //those displacements will be sent to the BEM code in the non-linear iterative solver
- // mapping nodeTag (of node on the boundary) into (u,v) displacements of the node
- retrieveBoundaryDisplacements(boundaryDisplacementMap, full_d_vector, nodeIndexMap, groups);
-
std::vector<std::string> names;
gmsh::model::list(names);
+ /*-------------POST PROCESSING---------------------*/
if(postProcessing)
{
- /*------------- ROTATION VALIDATION: APPLY 90 degree rotation-------------*/
- //elementData* element = &FEM_kinematicsMap[6];
- //std::cout << "hello: " << element->theta << "\n";
- //std::cout << "pl: " << element->pl << "\n";
- /*for(size_t i = 0; i < FEM_nodeTags.size(); i++)
- {
- std::vector<double> coord, parametricCoord;
- int dim, tag;
- double angle = M_PI /4;
- gmsh::model::mesh::getNode(FEM_nodeTags[i], coord, parametricCoord, dim, tag);
- std::vector<double> total_disp = {cos(angle)*(coord[0] + full_d_vector[2*i])-sin(angle)*(coord[1] + full_d_vector[2*i + 1]), cos(angle)*(coord[1] + full_d_vector[2*i + 1])+sin(angle)*(coord[0] + full_d_vector[2*i]), coord[2]};
- full_d_vector[2*i] = total_disp[0] - coord[0];
- full_d_vector[2*i+1] = total_disp[1] - coord[1];
- //std::cout << "nodetag: " << FEM_nodeTags[i] << ", new u: " << full_d_vector[2*i] << ", new v: " << full_d_vector[2*i+1] << "\n";
- }
- updateKinematics(FEM_kinematicsMap, full_d_vector, nodeIndexMap, FEM_elementTags);*/
-
- // retrieving maximal nodal displacement
+ // retrieving and displaying maximal nodal displacement
double max_disp = 0;
double max_u = 0;
double max_v = 0;
@@ -423,63 +364,19 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
std::cout << "maximal nodal displacement (u,v) = (" << max_u << "," << max_v << ") [m], at node " << max_disp_nodeTag << "\n";
/*--------------REACTION FORCES COMPUTATION------------------*/
- computeReactionForces(groups, nodeIndexMap, final_nodal_forces, keys);
+ computeReactionForces(groups, nodeIndexMap, final_nodal_forces, BCkeys);
/*--------------PLOTTING NODAL VALUES------------------*/
plotUxUy(FEM_nodeTags, full_d_vector, names);
/*-------------POST PROCESSING---------------*/
- /*-------------Computing and saving strains/stresses-----------------*/
-
-
- /*-------------------ELEMENTAL-NODAL VALUES-----------------------*/
- int viewTagEpsX = gmsh::view::add("nodal_eps_x");
- int viewTagEpsY = gmsh::view::add("nodal_eps_y");
- int viewTagEpsXY = gmsh::view::add("nodal_2eps_xy");
- int viewTagEpsZ = gmsh::view::add("nodal_eps_z");
- int viewTagSigX = gmsh::view::add("nodal_sig_x");
- int viewTagSigY = gmsh::view::add("nodal_sig_y");
- int viewTagSigXY = gmsh::view::add("nodal_sig_xy");
- int viewTagSigVM = gmsh::view::add("nodal_sig_VM");
-
- std::vector<std::vector<std::vector<double>>> dataEps(4);
- dataEps[0].resize(nbelems); // eps X
- dataEps[1].resize(nbelems); // eps Y
- dataEps[2].resize(nbelems); // eps XY
- dataEps[3].resize(nbelems); // eps Z
-
- std::vector<std::vector<std::vector<double>>> dataSig(4);
- dataSig[0].resize(nbelems); // Sig X
- dataSig[1].resize(nbelems); // Sig Y
- dataSig[2].resize(nbelems); // Sig XY
- dataSig[3].resize(nbelems); // Sig VM
- std::vector<std::size_t> elems2D(nbelems);
-
- // computing the elemental nodal data for each physical group of the domain
+ /*-------------Computing and saving elemental-nodal strains/stresses-----------------*/
double max_VM_stress = 0;
int max_VM_nodeTag;
- computeStressStrainElNodal(dim, tags, elems2D, H_matrix, nu, E, FEM_kinematicsMap, dataEps, dataSig, max_VM_stress, max_VM_nodeTag);
+ computeStressStrainElNodal(dim, tags, H_matrix, nu, E, FEM_kinematicsMap, nbelems, max_VM_stress, max_VM_nodeTag);
std::cout << "-----------------------------------------\n";
std::cout << "maximal nodal von Mises stress: " << max_VM_stress << " [Pa], at node " << max_VM_nodeTag << "\n";
-
- gmsh::view::addModelData(viewTagEpsX, 0, names[0], "ElementNodeData",
- elems2D, dataEps[0], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagEpsY, 0, names[0], "ElementNodeData",
- elems2D, dataEps[1], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagEpsXY, 0, names[0], "ElementNodeData",
- elems2D, dataEps[2], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagEpsZ, 0, names[0], "ElementNodeData",
- elems2D, dataEps[3], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigX, 0, names[0], "ElementNodeData",
- elems2D, dataSig[0], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigY, 0, names[0], "ElementNodeData",
- elems2D, dataSig[1], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigXY, 0, names[0], "ElementNodeData",
- elems2D, dataSig[2], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigVM, 0, names[0], "ElementNodeData",
- elems2D, dataSig[3], 1); // the last ,1 is important!
-
// retrieving work done by external forces
double externalWork = 0;
computeWorkDoneByExternalForces(externalWork, full_d_vector, final_nodal_forces);
@@ -497,7 +394,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
std::cout << "-----------------------------------------\n";
std::cout << "Total potential energy: " << strainEnergy - externalWork << " [J/m].\n";
}
- // representing vector displacement field and nodal forces
+ // representing vector displacement field and nodal forces as vectorial quantities
std::vector<std::vector<double>> dataU(FEM_nodeTags.size());
std::vector<std::vector<double>> dataF(FEM_nodeTags.size());
for (std::size_t i = 0; i < FEM_nodeTags.size(); ++i)
@@ -520,20 +417,29 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
if(postProcessing)
{
- // empeche que tout se plot l'un sur l'autre
- nbViews = 12 + nbViews; //hardcoded
+ nbViews = 12 + nbViews; // FEM code adds 12 views
+ // avoiding that all the views are superposed
for( int i = 0; i < nbViews; ++i){
std::string my_string = "View[" + std::to_string(i) + "].Visible";
gmsh::option::setNumber(my_string, 0);
}
- gmsh::option::setNumber("View[0].VectorType", 5); // plot automatiquement sous forme de déplacement
- gmsh::option::setNumber("View[0].RangeType", 3); // abscisses par step
+ gmsh::option::setNumber("View[0].VectorType", 5); // plot the displacement in deformed configuration
+ gmsh::option::setNumber("View[0].RangeType", 3);
}
double end_postPro = omp_get_wtime();
step_name = "post processing";
displayTime(start_postPro, end_postPro, step_name, display_time);
+ /*--------------RETRIEVING NODAL DISPLACEMENTS ON BEM-FEM BOUNDARY------------*/
+ // those displacements will be sent to the BEM code in the non-linear iterative solver
+ // mapping nodeTag (of node on the boundary) into (u,v) displacements of the node
+ retrieveBoundaryDisplacements(boundaryDisplacementMap, full_d_vector, nodeIndexMap, groups);
+
+ // if mesh untangler is activated, all nodal coordinates are updated according to their final
+ // nodal displacements. It is useful when coupled to the BEM solver in order to untangle
+ // the BEM domain and compute the potential field in the whole BEM domain adapted to the
+ // displacement of the FEM domain.
if(untangleMesh)
{
for(size_t i = 0; i < FEM_nodeTags.size(); i++)
@@ -546,10 +452,11 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
}
+ // it is only useful when using "large_rotation_validation.geo" file.
+ // printing the displacement of "point A" (see report for graph).
if(validation)
{
- // point A
- std::string groupname = "point_A"; // we consider the domain as containing the FEM 2d model
+ std::string groupname = "point_A";
int dim = groups[groupname].first;
int tag = groups[groupname].second;
std::vector<int> tags;
@@ -563,17 +470,13 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
}
}
-
+// solves the linear elastic problem in the FEM domain.
+// "electrostaticPressure" contains the pressure resulting from the electric field computation on the
+// BEM_FEM_boundary. "nbViews" is the current number of views already allocated in the gmsh window.
void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
{
- /*--------------INIT----------------*/
-
- int max_threads = omp_get_max_threads();
-
-
/*--------Integration rule used in the code----------*/
- // could be useful to pass it as an argument in the .geo file
- std::string integration_rule = "CompositeGauss4"; // tested with other methods, OK too.
+ std::string integration_rule = "CompositeGauss4";
/*--------get physical groups--------------*/
std::map<std::string, std::pair<int, int>> groups;
@@ -593,63 +496,37 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
int dim = groups[groupname].first;
int tag = groups[groupname].second;
- //loop over the entities of the domain to retrieve all nodeTags by taking care not to take some nodes twice
- //using sort and set_union methods
- std::vector<std::size_t> FEM_nodeTags;
- std::vector<double> FEM_nodecoord;
- std::vector<double> FEM_nodeparametricCoord;
- // Getting entities of the domain to fill FEM_nodeTags, en faire une fonction ?
+ // Getting entities of the domain to fill FEM_nodeTags
std::vector<int> tags;
gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
- std::vector<std::vector<std::size_t>> tmp_nodeTags(tags.size()); // will store the nodeTags associated to one given entity of the domain
- gmsh::model::mesh::getNodes(tmp_nodeTags[0], FEM_nodecoord, FEM_nodeparametricCoord, dim, tags[0], true);
- std::sort(tmp_nodeTags[0].begin(), tmp_nodeTags[0].end()); // need to sort the vectors in order to merge them afterwards
- for (std::size_t i = 1; i < tags.size(); i++)
- { // in the case the domain contains multiple entities
- gmsh::model::mesh::getNodes(tmp_nodeTags[i], FEM_nodecoord, FEM_nodeparametricCoord, dim, tags[i], true);
- std::sort(tmp_nodeTags[i].begin(), tmp_nodeTags[i].end());
- std::set_union(tmp_nodeTags[i-1].begin(), tmp_nodeTags[i-1].end(),
- tmp_nodeTags[i].begin(), tmp_nodeTags[i].end(), std::back_inserter(FEM_nodeTags));
-
- }
- if(tags.size() == 1){ // if the domain contains one single entity
- FEM_nodeTags = tmp_nodeTags[0];
- }
- else{ // little trick for more than 2 domains
- std::sort(FEM_nodeTags.begin(), FEM_nodeTags.end());
- auto last = std::unique(FEM_nodeTags.begin(), FEM_nodeTags.end());
- FEM_nodeTags.erase(last, FEM_nodeTags.end());
- }
-
+ std::vector<std::size_t> FEM_nodeTags;
+ std::map<int,std::pair<double,double>> FEM_nodeCoordsMap; // nodeTag -> (x,y) (useless in this linear case)
+ RetrieveFEMNodeTagsAndCoords(FEM_nodeTags, FEM_nodeCoordsMap, dim, tags);
/*-------map between nodeTag and index of first nodal displacement associated to the node (second one is the same +1)------*/
std::map<int, int> nodeIndexMap;
for (std::size_t i = 0; i < FEM_nodeTags.size(); ++i)
{
- int nodeTag = FEM_nodeTags[i];// use i à la place de nodeTag pour traiter cas ou le nodeTag n'est pas une suite continue
+ int nodeTag = FEM_nodeTags[i];
int node_first_index = 2*i;
nodeIndexMap[nodeTag] = node_first_index;
}
/*-----------get physical properties of domain----------------------*/
- std::vector<std::string> keys;
- gmsh::onelab::getNames(keys, "(Volumic Forces|Materials).+");
double E = 0; // Young modulus
double nu = 0; // Poisson ratio
double rho = 0; // volumic mass
double bx = 0; // volumic force along x axis
double by = 0; // volumic force along y axis
-
- getPhysicalProperties(keys, E, nu, rho, bx, by);
+ getPhysicalProperties(E, nu, rho, bx, by);
/*-----------------H matrix involved in stress computation------------------*/
Eigen::Matrix<double, 3, 3> H_matrix = computeHmatrix(E, nu);
- //double startbis = omp_get_wtime();
- //std::cout << "preliminary: " << startbis-start << "\n";
/*------------computing K matrix and f vector------------------------*/
// we first consider the K matrix and the volumic part of the f vector
+ // K matrix declared as sparse matrix for efficiency.
Eigen::SparseMatrix<double> full_K_matrix(2*FEM_nodeTags.size(),2*FEM_nodeTags.size());
std::vector<double> full_f_vector(2*FEM_nodeTags.size());
for (size_t i = 0; i < full_f_vector.size(); i++)
@@ -657,72 +534,21 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
full_f_vector[i] = 0;
}
- std::vector<std::vector<Eigen::Triplet <double>>> my_k_list(max_threads);
- std::vector<std::vector<std::pair<int,double>>> thread_f_vector(max_threads);
-
- // Looping over entities and compute full K (under tripletlist form) and F
- assembleKlistVolumicF(H_matrix, rho, bx, by, integration_rule, dim, tags, nodeIndexMap, my_k_list, thread_f_vector);
-
-
- std::vector<Eigen::Triplet <double>> k_tripletList = my_k_list[0];
-
- if(max_threads > 1)
- {
- for (int i = 1; i < max_threads; i++)
- {
- k_tripletList.insert(k_tripletList.end(), my_k_list[i].begin(), my_k_list[i].end());
- }
- }
- full_K_matrix.setFromTriplets(k_tripletList.begin(), k_tripletList.end());
- for(int i = 0; i < max_threads; i++){
- std::vector<std::pair<int,double>> &tmp_vector = thread_f_vector[i];
- for(auto &p : tmp_vector) {
- full_f_vector[p.first] += p.second;
- }
- }
-
- //double start2 = omp_get_wtime();
- //std::cout << "time for K assembly: " << start2-startbis << "\n";
-
+ // Looping over entities and compute full K and volumic F
+ assembleKVolumicF(H_matrix, rho, bx, by, integration_rule, dim, tags, nodeIndexMap, full_K_matrix, full_f_vector);
/*-----------------BOUNDARY CONDITIONS------------*/
/*-----------get group_names related to BC's----------------------*/
- gmsh::onelab::getNames(keys, "(Boundary Conditions).+");
+ std::vector<std::string> BCkeys;
+ gmsh::onelab::getNames(BCkeys, "(Boundary Conditions).+");
/*-----------first focus on neumann BC's (surface traction) into full_f_vector------------*/
- includeNeumannBC(keys, integration_rule, nodeIndexMap, groups, full_f_vector);
-
- //double start3 = omp_get_wtime();
- //std::cout << "time for neumann: " << start3 - start2 << "\n";
+ includeNeumannBC(BCkeys, integration_rule, nodeIndexMap, groups, full_f_vector);
/* ELECTROSTATIC PRESSURE */
- if(electrostaticPressure.size() != 0)
+ if(electrostaticPressure.size() != 0) // size is only greater than zero if there was already a BEM solver
{
- std::string groupnameBoundary = "BEM_FEM_boundary";
- int dimBoundary = groups[groupnameBoundary].first;
- int tagBoundary = groups[groupnameBoundary].second;
- std::vector<int> tagsBoundary;
- gmsh::model::getEntitiesForPhysicalGroup(dimBoundary, tagBoundary, tagsBoundary);
- std::vector<int> elementTypes;
- std::vector<std::vector<std::size_t>> elementTags;
- std::vector<std::vector<std::size_t>> elnodeTags;
-
- std::map<int,std::pair<double,double>> nodeCoordsMap;
- //filling the map with : nodeTag -> (x,y)
- for(std::size_t i=0; i< tags.size(); ++i)
- {
- std::vector<std::size_t> entity_nodeTags;
- std::vector<double> entity_nodecoord;
- std::vector<double> entity_nodeparametricCoord;
- gmsh::model::mesh::getNodes(entity_nodeTags, entity_nodecoord, entity_nodeparametricCoord, dim, tags[i], true);
- for ( std::size_t j=0; j < entity_nodeTags.size(); j++){
- nodeCoordsMap[entity_nodeTags[j]] = std::pair<double, double>(entity_nodecoord[3*j], entity_nodecoord[3*j+1]);
- }
- }
- applyElecPressure(integration_rule, nodeCoordsMap, nodeIndexMap, groups, electrostaticPressure, full_f_vector);
+ applyElecPressure(integration_rule, FEM_nodeCoordsMap, nodeIndexMap, groups, electrostaticPressure, full_f_vector);
}
-
- //double start4 = omp_get_wtime();
- //std::cout << "time for elec pressure: " << start4 - start3 << "\n";
/*------------FOCUS ON DIRICHLET BC'S---------------*/
// new_DOFindices vector will contain the new indices of the K matrix after the rows/columns corresponding
@@ -736,29 +562,15 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
std::vector<double> dir_BC(2*FEM_nodeTags.size());
// Fill the two vectors defined just above
- fillDOFindicesDirBC(keys, nodeIndexMap, groups, new_DOFindices, dir_BC);
-
int number_removed_lines = 0;
- for (std::size_t j = 0; j < new_DOFindices.size(); ++j){
- if(new_DOFindices[j] < 0){ // line has to be removed
- number_removed_lines++;
- }
- else{
- new_DOFindices[j] -= number_removed_lines;
- }
- }
-
- //double start3 = omp_get_wtime();
- //std::cout << "time for boundary conditions: " << start3 - start2 << "\n";
+ fillDOFindicesDirBC(BCkeys, nodeIndexMap, groups, new_DOFindices, dir_BC, number_removed_lines);
/*------------FILLING THE REDUCED K MATRIX AND F VECTOR--------------*/
- //the correction term on the f vector comes from the "FEM" theoretical course (JP Ponthot), Chap2 Slide 38/50
+ //the correction term on the f vector comes from the "Finite Element Method" theoretical course (JP Ponthot, ULiege), Chap2 Slide 38/50
Eigen::SparseMatrix<double> reduced_K_matrix(2*FEM_nodeTags.size() - number_removed_lines,2*FEM_nodeTags.size() - number_removed_lines);
std::vector<double> reduced_f_vector(2*FEM_nodeTags.size() - number_removed_lines);
- std::list<Eigen::Triplet <double>> reduced_k_tripletList;
-
- fillReducedKf(full_f_vector, full_K_matrix, new_DOFindices, dir_BC, reduced_K_matrix, reduced_f_vector, reduced_k_tripletList);
+ fillReducedKf(full_f_vector, full_K_matrix, new_DOFindices, dir_BC, reduced_K_matrix, reduced_f_vector);
/*--------------SOLVING THE REDUCED SYSTEM--------------------*/
Eigen::VectorXd final_reduced_d(2*FEM_nodeTags.size() - number_removed_lines);
@@ -769,15 +581,11 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
}
// solve Kd=f
- Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
+ Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
solver.compute(reduced_K_matrix);
final_reduced_d = solver.solve(final_reduced_f);
- //double start4 = omp_get_wtime();
- //std::cout << "time to solve system: " << start4 -start3 << "\n";
-
-
/*-----------RECONSTRUCTING FULL NODAL VALUES--------------*/
std::vector<double> full_d_vector(2*FEM_nodeTags.size());
for (std::size_t j = 0; j < full_d_vector.size(); ++j){
@@ -803,12 +611,6 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
nodal_forces[i] = nodalForces(i);
}
- /*--------------REACTION FORCES COMPUTATION------------------*/
- computeReactionForces(groups, nodeIndexMap, nodal_forces, keys);
-
- //double start5 = omp_get_wtime();
- //std::cout << "time for rest: " << start5 - start4 << "\n";
-
/*--------------PLOTTING AND SAVING NODAL VALUES------------------*/
std::vector<std::string> names;
gmsh::model::list(names);
@@ -836,49 +638,38 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
}
}
+ std::cout << "The FEM domain contains " << nbelems << " element(s) and " << FEM_nodeTags.size() << " nodes\n";
+
+ // retrieving and displaying maximal nodal displacement
+ double max_disp = 0;
+ double max_u = 0;
+ double max_v = 0;
+ int max_disp_nodeTag = 0;
+
+ for(size_t i = 0; i < FEM_nodeTags.size(); i++)
+ {
+ double tmp_index = nodeIndexMap[FEM_nodeTags[i]];
+ double tmp_disp = full_d_vector[tmp_index]*full_d_vector[tmp_index] + full_d_vector[tmp_index + 1]*full_d_vector[tmp_index + 1];
+ if (tmp_disp > max_disp)
+ {
+ max_disp = tmp_disp;
+ max_disp_nodeTag = FEM_nodeTags[i];
+ max_u = full_d_vector[tmp_index];
+ max_v = full_d_vector[tmp_index + 1];
+ }
+ }
+ std::cout << "-----------------------------------------\n";
+ std::cout << "maximal nodal displacement (u,v) = (" << max_u << "," << max_v << ") [m], at node " << max_disp_nodeTag << "\n";
+
+ /*--------------REACTION FORCES COMPUTATION------------------*/
+ computeReactionForces(groups, nodeIndexMap, nodal_forces, BCkeys);
/*-------------------ELEMENTAL-NODAL VALUES-----------------------*/
- int viewTagEpsX = gmsh::view::add("nodal_eps_x");
- int viewTagEpsY = gmsh::view::add("nodal_eps_y");
- int viewTagEpsXY = gmsh::view::add("nodal_2eps_xy");
- int viewTagEpsZ = gmsh::view::add("nodal_eps_z");
- int viewTagSigX = gmsh::view::add("nodal_sig_x");
- int viewTagSigY = gmsh::view::add("nodal_sig_y");
- int viewTagSigXY = gmsh::view::add("nodal_sig_xy");
- int viewTagSigVM = gmsh::view::add("nodal_sig_VM");
-
- std::vector<std::vector<std::vector<double>>> dataEps(4);
- dataEps[0].resize(nbelems); // eps X
- dataEps[1].resize(nbelems); // eps Y
- dataEps[2].resize(nbelems); // eps XY
- dataEps[3].resize(nbelems); // eps Z
-
- std::vector<std::vector<std::vector<double>>> dataSig(4);
- dataSig[0].resize(nbelems); // Sig X
- dataSig[1].resize(nbelems); // Sig Y
- dataSig[2].resize(nbelems); // Sig XY
- dataSig[3].resize(nbelems); // Sig VM
- std::vector<std::size_t> elems2D(nbelems);
-
- // computing the elemental nodal data for each physical group of the domain
- computeStressStrainLinearWay(dim, tags, full_d_vector, elems2D, H_matrix, nu, E, nodeIndexMap, dataEps, dataSig);
-
- gmsh::view::addModelData(viewTagEpsX, 0, names[0], "ElementNodeData",
- elems2D, dataEps[0], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagEpsY, 0, names[0], "ElementNodeData",
- elems2D, dataEps[1], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagEpsXY, 0, names[0], "ElementNodeData",
- elems2D, dataEps[2], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagEpsZ, 0, names[0], "ElementNodeData",
- elems2D, dataEps[3], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigX, 0, names[0], "ElementNodeData",
- elems2D, dataSig[0], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigY, 0, names[0], "ElementNodeData",
- elems2D, dataSig[1], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigXY, 0, names[0], "ElementNodeData",
- elems2D, dataSig[2], 1); // the last ,1 is important!
- gmsh::view::addModelData(viewTagSigVM, 0, names[0], "ElementNodeData",
- elems2D, dataSig[3], 1); // the last ,1 is important!
+ double max_VM_stress = 0;
+ int max_VM_nodeTag;
+ computeStressStrainLinearWay(dim, tags, full_d_vector, H_matrix, nu, E, nodeIndexMap, nbelems, max_VM_stress, max_VM_nodeTag);
+ std::cout << "-----------------------------------------\n";
+ std::cout << "maximal nodal von Mises stress: " << max_VM_stress << " [Pa], at node " << max_VM_nodeTag << "\n";
// representing vector displacement field and nodal forces
int viewTagU = gmsh::view::add("u");
@@ -904,12 +695,9 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews)
gmsh::view::addModelData(viewTagF, 0, names[0], "NodeData",
FEM_nodeTags, dataF); //independent of time here
-
-
- // empeche que tout se plot l'un sur l'autre
- int nb_views = 12 + nbViews; //hardcoded
+ int nb_views = 12 + nbViews; //FEM code adds 12 views
for( int i = 0; i < nb_views; ++i){
- std::string my_string = "View[" + std::to_string(i) + "].Visible";
+ std::string my_string = "View[" + std::to_string(i) + "].Visible"; // avoids that all views are represented on top of each other
gmsh::option::setNumber(my_string, 0);
}
--
GitLab