Skip to content
Snippets Groups Projects
Commit f04e4330 authored by Denis Louis's avatar Denis Louis
Browse files

begin kinematics implementation

parent c77b0362
No related branches found
No related tags found
1 merge request!8SECOND PROGRESS DEADLINE
Hide whitespace changes
Inline Side-by-side
srcs/FEM/functionsFEM.cpp
+ 106
− 14
View file @ f04e4330
......@@ -4,6 +4,8 @@
#include <gmsh.h> // gmsh main header
#endif
#include "functionsFEM.hpp"
#include <iostream> // for std::cout
#include <map>
#include <sstream> // for std::stringstream
......@@ -718,7 +720,7 @@ void Compute_Stress_Strain_Elemental_Nodal(const std::vector<int> & elementTypes
}
}
}
void applyElecPressure(const std::string & integration_rule, std::map<int, int> & node_map, std::map<std::string, std::pair<int, int>> & groups, std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector){
void applyElecPressure(const std::string & integration_rule, std::map<int,std::pair<double,double>> & nodeCoordsMap, std::map<int, int> & node_map, 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;
......@@ -729,19 +731,6 @@ void applyElecPressure(const std::string & integration_rule, std::map<int, int>
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]);
}
}
for(std::size_t i=0; i< tags.size(); ++i)
{
gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, dim, tags[i]);
......@@ -814,4 +803,107 @@ void applyElecPressure(const std::string & integration_rule, std::map<int, int>
}
}
}
}
int initKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags)
{
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
std::vector<std::vector<std::size_t>> elnodeTags;
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();
// getting element properties, in particular: number of nodes
std::string elementName;
int elDim, order, numNodes, numPrimaryNodes;
std::vector<double> localNodeCoord;
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;
FEM_kinematicsMap[elementTags[j][k]] = element;
}
}
}
return nbelems;
}
void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector, std::map<int, int> & node_map, int & dim, std::vector<int> & tags)
{
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
std::vector<std::vector<std::size_t>> elnodeTags;
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)
{
for(std::size_t k = 0; k < elementTags[j].size(); k++)
{
elementData* element = & FEM_kinematicsMap[elementTags[j][k]];
for(size_t l = 0; l < element->nodes.size(); l++)
{
size_t tmp_nodeTag = element->nodes[l];
int node_index = node_map[tmp_nodeTag];
element->u[l] = full_d_vector[node_index];
element->v[l] = full_d_vector[node_index + 1];
}
// not finished yet !
}
}
}
}
\ No newline at end of file
srcs/FEM/functionsFEM.hpp
+ 27
− 1
View file @ f04e4330
......@@ -5,6 +5,30 @@
#include <list>
#include <Eigen/Sparse> // Eigen library for sparse matrices
struct elementData{
// will store all the kinematics data
std::vector<size_t> nodes; // will contain the tags of the nodes
std::vector<double> x; // will contain the x global coordinate of the nodes
std::vector<double> y; // will contain the y global coordinate of the nodes
double xc; // x global coordinate of center
double yc; // y global coordinate of center
std::vector<double> xi; // x local coordinate of nodes
std::vector<double> yi; // y local coordinate of nodes
std::vector<double> u; // x global displacement of the nodes
std::vector<double> v; // y global displacement of the nodes
double uc; // x global displacement of center
double vc; // y global displacement of center
double theta; // angle of rotation of element
std::vector<double> ui; // x local displacement of nodes
std::vector<double> vi; // y local displacement of nodes
};
Eigen::Matrix<double, 3, 3> Compute_H_matrix(const double E, const double nu);
void Assembly_K_triplet(const int i, const int el, const std::vector<std::vector<std::size_t>> & nodeTags, const int numNodes, std::map<int, int> & node_map, std::list<Eigen::Triplet <double>> & k_tripletList, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix, std::vector<double> & f_vector);
void Assembly_f_thread(const int i, const int el, const std::vector<std::vector<std::size_t>> & nodeTags, const int numNodes, std::map<int, int> & node_map, std::list<std::pair<int,double>> & thread_f_vector, std::vector<double> & f_vector);
......@@ -18,6 +42,8 @@ void Fill_K_F_Reduced(const std::vector<double> & full_f_vector, const Eigen::Sp
void Reaction_Forces(double & Rx, double & Ry, std::map<std::string, std::pair<int, int>> & groups, std::map<int, int> & node_map, Eigen::VectorXd & nodal_forces, const std::vector<std::string> keys);
std::vector<std::string> Plot_Save_ux_uy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector, const int nbViews);
void Compute_Stress_Strain_Elemental_Nodal(const std::vector<int> & elementTypes, const std::vector<std::vector<std::size_t>> &elementTags, const std::vector<std::vector<std::size_t>> &elnodeTags, 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> & node_map, std::vector<std::vector<double>> & data3, std::vector<std::vector<double>> & data4, std::vector<std::vector<double>> & data5, std::vector<std::vector<double>> & data6, std::vector<std::vector<double>> & data7, std::vector<std::vector<double>> & data8, std::vector<std::vector<double>> & data9, std::vector<std::vector<double>> & data10, int & k);
void applyElecPressure(const std::string & integration_rule, std::map<int, int> & node_map, std::map<std::string, std::pair<int, int>> & groups, std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector);
void applyElecPressure(const std::string & integration_rule, std::map<int,std::pair<double,double>> & nodeCoordsMap, std::map<int, int> & node_map, std::map<std::string, std::pair<int, int>> & groups, std::map<int, double> &electrostaticPressure, std::vector<double> & full_f_vector);
int initKinematics(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> & node_map, int & dim, std::vector<int> & tags);
int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressure, const int nbViews);
\ No newline at end of file
srcs/FEM/solverFEM.cpp
+ 55
− 23
View file @ f04e4330
......@@ -56,19 +56,28 @@ int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
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
//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));
......@@ -78,6 +87,27 @@ int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
FEM_nodeTags = tmp_nodeTags[0];
}
/*-----get elements related to FEM--------------*/
std::map<int,elementData> FEM_kinematicsMap; // map between elementTag and all the kinematics associated to this element
// map is initialized, see functionsFEM.hpp for structure of elementData type
int nbelems = initKinematics(FEM_kinematicsMap, FEM_nodeCoordsMap, dim, tags);
elementData myTest = FEM_kinematicsMap[533]; // 533 hardcoded, simply choose one elementTag of the geometry
std::cout << "my test nodes: ";
for(size_t i = 0; i < myTest.nodes.size(); i++)
{
std::cout << myTest.nodes[i] << ", ";
}
std::cout << "\n";
std::cout << "my test xi: ";
for(size_t i = 0; i < myTest.xi.size(); i++)
{
std::cout << myTest.xi[i] << ", ";
}
std::cout << "\n";
std::cout << "my test yc: " << myTest.yc << "\n";
/*-------map between nodeTag and index of first nodal displacement associated to the node (second one is the same +1)------*/
......@@ -158,7 +188,7 @@ int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
/* ELECTROSTATIC PRESSURE */
if(electrostaticPressure.size() != 0)
applyElecPressure(integration_rule, node_map, groups, electrostaticPressure, full_f_vector);
applyElecPressure(integration_rule, FEM_nodeCoordsMap, node_map, groups, electrostaticPressure, full_f_vector);
//double start4 = omp_get_wtime();
//std::cout << "time for elec pressure: " << start4 - start3 << "\n";
......@@ -227,6 +257,25 @@ int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
}
}
/*--------- UPDATING KINEMATICS---------------*/
updateKinematics(FEM_kinematicsMap, full_d_vector, node_map, dim, tags);
myTest = FEM_kinematicsMap[24]; // 533 hardcoded, simply choose one elementTag of the geometry
std::cout << "my test nodes: ";
for(size_t i = 0; i < myTest.nodes.size(); i++)
{
std::cout << myTest.nodes[i] << ", ";
}
std::cout << "\n";
std::cout << "my test ui: ";
for(size_t i = 0; i < myTest.u.size(); i++)
{
std::cout << myTest.u[i] << " " << full_d_vector[node_map[myTest.nodes[i]]] << ", ";
}
std::cout << "\n";
std::cout << "my test vc: " << myTest.vc << "\n";
/*--------------COMPUTING NODAL FORCES-------------*/
Eigen::VectorXd eigen_d(2*FEM_nodeTags.size());
for(size_t i = 0 ; i < 2*FEM_nodeTags.size(); ++i){
......@@ -251,27 +300,6 @@ int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
/*-------------POST PROCESSING---------------*/
/*-------------Computing and saving strains/stresses-----------------*/
// get all the elements in the domain
groupname = "domain"; // we consider the domain as containing the FEM 2d model
dim = groups[groupname].first;
tag = groups[groupname].second;
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;
// compute total number of elements in the domain
int nbelems=0;
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();
}
}
/*-------------------ELEMENTAL-NODAL VALUES-----------------------*/
// éventuellement faire un tableau ici qui regroupe tout le bazar pour ne pas avoir la même ligne 7 fois ...
......@@ -292,6 +320,10 @@ int solverFEM(int argc, char **argv, std::map<int, double> &electrostaticPressur
std::vector<std::vector<double>> data9(nbelems);
std::vector<std::vector<double>> data10(nbelems);
std::vector<std::size_t> elems2D(nbelems);
std::vector<int> elementTypes;
std::vector<std::vector<std::size_t>> elementTags;
std::vector<std::vector<std::size_t>> elnodeTags;
// computing the elemental nodal data for each physical group of the domain
int k = 0; // current index of data vectors, k is incremented in Compute_Stress_Strain_Elemental_Nodal at each element
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment