G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsPeriodicOverlap.hpp
Go to the documentation of this file.
1
15namespace gismo
16{
17
18template <class T>
20{
21 calculate(this->m_options.getInt("parametrizationMethod"));
22}
23
24template<class T>
25void gsPeriodicOverlap<T>::calculate(const size_t paraMethod)
26{
27 size_t n = this->m_mesh.getNumberOfInnerVertices();
28 size_t N = this->m_mesh.getNumberOfVertices();
29
30 Neighbourhood neighbourhood(this->m_mesh, paraMethod);
31
32 this->initParameterPoints();
33
34 // Construct the twins.
35 constructTwins();
36
37 // Solve.
38 constructAndSolveEquationSystem(neighbourhood, n, N);
39}
40
41
42
43
44template <class T>
45void gsPeriodicOverlap<T>::constructTwinsBetween(size_t& currentNrAllVertices,
46 std::list<size_t> vertexIndices,
47 size_t from, size_t to,
48 bool rightHandSide)
49{
50 // TODO: The whiles do not check if the sought member is indeed in
51 // the list (danger of an infinite loop).
52
53 // Rotate the vertexIndices so as to start from from.
54 while(vertexIndices.front() != from)
55 {
56 vertexIndices.push_back(vertexIndices.front());
57 vertexIndices.pop_front();
58 }
59
60 // Push the corresponding pairs to the twin vector.
61 // Note that if we would do std::prev on begin(), we should not dereference.
62 // The error is visible only with -DGISMO_WITH_XDEBUG=ON.
63 for(std::list<size_t>::const_iterator it=vertexIndices.begin();
64 it == vertexIndices.begin() || *std::prev(it) != to;
65 ++it)
66 {
67 size_t twin = findTwin(*it);
68 if(rightHandSide)
69 m_twins.push_back(std::pair<size_t, size_t>(twin, ++currentNrAllVertices));
70 else
71 m_twins.push_back(std::pair<size_t, size_t>(++currentNrAllVertices, twin));
72 }
73
74}
75
76template<class T>
78{
79 // vertex with parameter v = 0 and lowest u value
80 gsVertexHandle uMinv0 = this->m_mesh.getVertex(this->m_indicesV0.front());
81
82 // vertex with parameter v = 0 and highest u value
83 gsVertexHandle uMaxv0 = this->m_mesh.getVertex(this->m_indicesV0.back());
84
85 // vertex with parameter v = 1 and lowest u value
86 gsVertexHandle uMinv1 = this->m_mesh.getVertex(this->m_indicesV1.front());
87
88 // vertex with parameter v = 1 and highest u value
89 gsVertexHandle uMaxv1 = this->m_mesh.getVertex(this->m_indicesV1.back());
90
91 const std::list<size_t> vertexIndices = m_overlapHEM.getBoundaryVertexIndices();
92 size_t currentNrAllVertices = this->m_mesh.getNumberOfVertices();
93
94 // Construct the twins on the right boundary of the overlap.
95 constructTwinsBetween(currentNrAllVertices, vertexIndices, uMinv1, uMinv0, true);
96
97 // Construct the twins on the left boundary of the overlap.
98 constructTwinsBetween(currentNrAllVertices, vertexIndices, uMaxv0, uMaxv1, false);
99}
100
101template <class T>
103 const size_t n,
104 const size_t N)
105{
106 size_t numTwins = m_twins.size();
107 gsMatrix<T> LHS(N + numTwins, N + numTwins);
108 gsMatrix<T> RHS(N + numTwins, 2);
109 // prevent Valgrind warnings
110 LHS.setZero();
111 RHS.setZero();
112 std::vector<T> lambdas;
113
114 // interior points
115 for (size_t i = 0; i < n; i++)
116 {
117 lambdas = neighbourhood.getLambdas(i);
118 updateLambdasWithTwins(lambdas, i+1);
119
120 for (size_t j = 0; j < N + numTwins; j++)
121 LHS(i, j) = ( i==j ? (T)(1) : -lambdas[j] );
122 }
123
124 // points on the lower and upper boundary
125 for (size_t i=n; i<N; i++)
126 {
127 LHS(i, i) = (T)(1);
128 RHS.row(i) = this->m_parameterPoints[i];
129 }
130
131 // points on the overlap
132 for (size_t i=N; i<N+numTwins; i++)
133 {
134 size_t first = m_twins[i-N].first-1;
135 size_t second = m_twins[i-N].second-1;
136
137 LHS(i, first) = (T)( 1);
138 LHS(i, second) = (T)(-1);
139
140 RHS(i, 0) = (T)(-1);
141 RHS(i, 1) = (T)( 0);
142 }
143
144 gsEigen::PartialPivLU<typename gsMatrix<T>::Base> LU = LHS.partialPivLu();
145 gsMatrix<T> sol = LU.solve(RHS);
146 for (size_t i = 0; i < N; i++)
147 {
148 this->m_parameterPoints[i] << sol(i, 0), sol(i, 1);
149 }
150}
151
152template <class T>
154 size_t vertexId) const
155{
156 lambdas.reserve(lambdas.size() + m_twins.size());
157 for(size_t i=0; i<m_twins.size(); i++)
158 lambdas.push_back(0);
159
160 // Determine, whether vertexId is on the left or right side of the overlap.
161 bool isLeft = false;
162 bool isRight = false;
163
164 for(auto it=m_twins.begin(); it!=m_twins.end(); ++it)
165 {
166 if(it->first == vertexId)
167 {
168 isRight = true;
169 break;
170 }
171 else if(it->second == vertexId)
172 {
173 isLeft = true;
174 break;
175 }
176 }
177
178
179 for(size_t i=0; i<m_twins.size(); i++)
180 {
181 size_t first=m_twins[i].first-1;
182 size_t second=m_twins[i].second-1;
183
184 // Left vertex swaps all its right neighbours.
185 if(isRight && first > second && lambdas[second] != 0)
186 {
187 lambdas[first] = lambdas[second];
188 lambdas[second] = 0;
189 }
190 // Right vertex swaps all its left neighbours
191 else if(isLeft && first < second && lambdas[first] != 0)
192 {
193 lambdas[second] = lambdas[first];
194 lambdas[first] = 0;
195 }
196 // Nothing happens to vertices that are not on the overlap.
197 }
198}
199
200template<class T>
202{
203 // Remember the vertices on the overlap boundaries.
204 std::vector<size_t> left, right;
205 for(auto it=m_twins.begin(); it!=m_twins.end(); ++it)
206 {
207 if(it->first < it->second)
208 left.push_back(it->first);
209 else
210 right.push_back(it->second);
211 }
212
213 typename gsPeriodicParametrization<T>::FlatMesh display(createExtendedFlatMesh(left, right));
214 return display.createRestrictedFlatMesh();
215}
216
217template<class T>
219 const std::vector<size_t>& left) const
220{
221 typedef typename gsParametrization<T>::Point2D Point2D;
222
223 gsMesh<T> midMesh;
224 midMesh.reserve(3 * this->m_mesh.getNumberOfTriangles(), this->m_mesh.getNumberOfTriangles(), 0);
225
226 for (size_t i = 0; i < this->m_mesh.getNumberOfTriangles(); i++)
227 {
228 size_t vInd[3];
229
230 // the indices of the current triangle vertices that are among left or right, respectively
231 std::vector<size_t> lVert, rVert;
232 for (size_t j = 1; j <= 3; ++j)
233 {
234 vInd[j-1] = this->m_mesh.getGlobalVertexIndex(j, i);
235 if(std::find(right.begin(), right.end(), vInd[j-1]) != right.end())
236 rVert.push_back(j-1);
237 if(std::find(left.begin(), left.end(), vInd[j-1]) != left.end())
238 lVert.push_back(j-1);
239 }
240
241 // Is the triangle inside overlap?
242 if(lVert.size() > 0 && rVert.size() > 0 && lVert.size() + rVert.size() == 3)
243 {
244 // Make two shifted copies of the triangle.
245 typename gsMesh<T>::VertexHandle mvLft[3], mvRgt[3];
246
247 for (size_t j=0; j<3; ++j)
248 {
249 const Point2D vertex = gsParametrization<T>::getParameterPoint(vInd[j]);
250 if(std::find(rVert.begin(), rVert.end(), j) != rVert.end())
251 {
252 mvLft[j] = midMesh.addVertex(vertex[0], vertex[1]);
253 mvRgt[j] = midMesh.addVertex(vertex[0]+(T)(1), vertex[1]);
254 }
255 else
256 {
257 mvLft[j] = midMesh.addVertex(vertex[0]-(T)(1), vertex[1]);
258 mvRgt[j] = midMesh.addVertex(vertex[0], vertex[1]);
259 }
260 }
261 midMesh.addFace(mvLft[0], mvLft[1], mvLft[2]);
262 midMesh.addFace(mvRgt[0], mvRgt[1], mvRgt[2]);
263 }
264 else
265 {
266 // Make just one triangle.
267 typename gsMesh<T>::VertexHandle mv[3];
268 for (size_t j=0; j<3; ++j)
269 {
270 const Point2D vertex = gsParametrization<T>::getParameterPoint(vInd[j]);
271 mv[j] = midMesh.addVertex(vertex[0], vertex[1]);
272 }
273 midMesh.addFace(mv[0], mv[1], mv[2]);
274 }
275 }
276 return midMesh.cleanMesh();
277}
278
279} // namespace gismo
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
Class Representing a triangle mesh with 3D vertices.
Definition gsMesh.h:32
gsMesh & cleanMesh()
reorders the vertices of all faces of an .stl mesh, such that only 1 vertex is used instead of #(adja...
Definition gsMesh.hpp:250
Class that maintains neighbourhood information of triangle mesh. Represents the neighbourhood propert...
Definition gsParametrization.h:286
const std::vector< T > & getLambdas(const size_t i) const
Get vector of lambdas.
Definition gsParametrization.hpp:282
const Point2D & getParameterPoint(size_t vertexIndex) const
Get parameter point Returns the parameter point with given vertex index.
Definition gsParametrization.hpp:160
gsMesh< T > createFlatMesh() const
Definition gsPeriodicOverlap.hpp:201
void constructAndSolveEquationSystem(const Neighbourhood &neighbourhood, const size_t n, const size_t N)
Analogous to the overloaded function from the parent class.
Definition gsPeriodicOverlap.hpp:102
void updateLambdasWithTwins(std::vector< T > &lambdas, size_t vertexId) const
Definition gsPeriodicOverlap.hpp:153
void compute()
Computes the periodic parametrization.
Definition gsPeriodicOverlap.hpp:19
gsMesh< T > createExtendedFlatMesh(const std::vector< size_t > &right, const std::vector< size_t > &left) const
Definition gsPeriodicOverlap.hpp:218
void constructTwins()
Construct the twins.
Definition gsPeriodicOverlap.hpp:77
void constructTwinsBetween(size_t &currentNrAllVertices, std::list< size_t > vertexIndices, size_t from, size_t to, bool rightHandSide)
Definition gsPeriodicOverlap.hpp:45
void calculate(const size_t paraMethod)
Definition gsPeriodicOverlap.hpp:25
Nested class for plotting flat meshes restricted to [0, 1]^2.
Definition gsPeriodicParametrization.h:33
gsMesh< T > createRestrictedFlatMesh() const
Trims the mesh to [0, 1]^2.
Definition gsPeriodicParametrization.hpp:154
A Point in T^d, with an index number.
Definition gsPoint.h:27
gsVertex class that represents a 3D vertex for a gsMesh.
Definition gsVertex.h:27
The G+Smo namespace, containing all definitions for the library.