Newer
Older
#ifndef COMPUTECENTER_H
#define COMPUTECENTER_H
#include "MeshFunctionsDefine.h"
#include "MeshApply.h"
#include "../MeshDataContainer/MeshDataContainer.h"
#include <array>
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
template <unsigned int dim, unsigned int Dimension, ComputationMethod Method = ComputationMethod::DEFAULT>
struct _ComputeCenters {
template <typename IndexType, typename Real, unsigned int ...Reserve>
static void compute(
MakeMeshDataContainer_t<Vertex<Dimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, Dimension>>& centers,
MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
auto& elemCenters = centers.template getDataByDim<dim>();
auto& subElemCenters = centers.template getDataByDim<dim - 1>();
for (IndexType i = 0; i < mesh.template getElements<dim>().size(); i++) {
auto& element = mesh.template getElements<dim>().at(i);
Real subElemCnt = 0;
for(auto& sub : element.getSubelements()){
elemCenters.at(i) += subElemCenters.at(sub.index);
subElemCnt++;
}
elemCenters.at(i) /= subElemCnt;
}
_ComputeCenters<dim + 1, Dimension, Method>::compute(centers, mesh);
}
};
template <unsigned int Dimension, ComputationMethod Method>
struct _ComputeCenters<Dimension, Dimension, Method>{
template <typename IndexType, typename Real, unsigned int ...Reserve>
static void compute(MakeMeshDataContainer_t<Vertex<Dimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, Dimension>>& centers,MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
auto& elemCenters = centers.template getDataByDim<Dimension>();
auto& subElemCenters = centers.template getDataByDim<Dimension - 1>();
for (IndexType i = 0; i < mesh.template getElements<Dimension>().size(); i++) {
auto& element = mesh.template getElements<Dimension>().at(i);
Real subElemCnt = 0;
IndexType tmpFaceIndex = element.getBoundaryElementIndex();
do {
elemCenters.at(i) += subElemCenters.at(tmpFaceIndex);
subElemCnt++;
tmpFaceIndex = mesh.getFaces()[tmpFaceIndex].getNextBElem(i);
} while (tmpFaceIndex != element.getBoundaryElementIndex());
elemCenters.at(i) /= subElemCnt;
}
}
};
template <unsigned int Dimension, ComputationMethod Method>
struct _ComputeCenters<1, Dimension, Method>{
template <typename IndexType, typename Real, unsigned int ...Reserve>
static void compute(
MakeMeshDataContainer_t<Vertex<Dimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, Dimension>>& centers,
MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
std::vector<Vertex<Dimension, Real>>& edgeCenters = centers.template getDataByDim<1>();
for (auto& edge : mesh.template getElements<1>()) {
edgeCenters.at(edge.getIndex()) = (mesh.template getElements<0>().at(edge.getVertexAIndex()) +
mesh.template getElements<0>().at(edge.getVertexBIndex())) * 0.5;
}
_ComputeCenters<2, Dimension, Method>::compute(centers, mesh);
}
};
/**
* @brief The _ComputeCenters<_Tp1, _Tp2, CalculationMethod::MESH_TESSELLATED> struct
* specialization of computation of centers for mesh which needs tessellation of faces
*/
template <>
struct _ComputeCenters<2, 3, ComputationMethod::TESSELLATED> {
template <typename IndexType, typename Real, unsigned int ...Reserve>
static void compute(
MakeMeshDataContainer_t<Vertex<3, Real>, make_custom_integer_sequence_t<unsigned int, 1, 3>>& centers,
MeshElements<3, IndexType, Real, Reserve...>& mesh){
DBGMSG("centers computation tessellated");
auto& elemCenters = centers.template getDataByDim<2>();
auto& subElemCenters = centers.template getDataByDim<2 - 1>();
for (IndexType i = 0; i < mesh.template getElements<2>().size(); i++) {
auto& element = mesh.template getElements<2>().at(i);
Real subElemCnt = 0;
for(auto& sub : element.getSubelements()){
elemCenters.at(i) += subElemCenters.at(sub.index);
subElemCnt++;
}
elemCenters.at(i) /= subElemCnt;
}
// calculate volume
for (IndexType i = 0; i < mesh.template getElements<2>().size(); i++) {
auto& element = mesh.template getElements<2>().at(i);
Vertex<3, Real> tempVert = {};
Real surfTotal = 0.0;
for(auto& sub : element.getSubelements()){
IndexType AI = mesh.getEdges().at(sub.index).getVertexAIndex();
IndexType BI = mesh.getEdges().at(sub.index).getVertexBIndex();
std::array<Vertex<3, Real>, 2> v = {elemCenters.at(i) - mesh.getVertices().at(AI), elemCenters.at(i) - mesh.getVertices().at(BI)};
std::array<Real, 2> norms;
grammSchmidt<2, 3, IndexType, Real>(v, norms);
Real surf = norms.at(0) * 0.5 * norms.at(1);
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
tempVert += subElemCenters.at(sub.index) * (surf * (2.0 / 3.0));
surfTotal += surf;
}
elemCenters.at(i) = (elemCenters.at(i) / 3.0) + (tempVert / surfTotal);
}
_ComputeCenters<2 + 1, 3, ComputationMethod::TESSELLATED>::compute(centers, mesh);
}
};
template <ComputationMethod Method, unsigned int Dimension,typename IndexType, typename Real, unsigned int ...Reserve>
MakeMeshDataContainer_t<Vertex<Dimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, Dimension>>
ComputeCenters(MeshElements<Dimension, IndexType, Real, Reserve...>& mesh){
MakeMeshDataContainer_t<Vertex<Dimension, Real>, make_custom_integer_sequence_t<unsigned int, 1, Dimension>> centers(mesh);
_ComputeCenters<1, Dimension, Method>::compute(centers, mesh);
return centers;
}
#endif // COMPUTECENTER_H