G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsC1SurfGluingDataVisitor.h
Go to the documentation of this file.
1
14#pragma once
15
16namespace gismo
17{
18
19
20template <class T>
21class gsC1SurfGluingDataVisitor
22{
23public:
24
25 gsC1SurfGluingDataVisitor()
26 {
27 }
28
29 void initialize(const gsBasis<T> & basis,
30 gsQuadRule<T> & rule)
31 {
32 gsVector<index_t> numQuadNodes( basis.dim() );
33 for (index_t i = 0; i < basis.dim(); ++i) // to do: improve
34 numQuadNodes[i] = 2 * basis.degree(i) + 1;
35
36 // Setup Quadrature
37 rule = gsGaussRule<T>(numQuadNodes);// NB!
38 }
39
40 // Evaluate on element.
41 inline void evaluate(gsMatrix<T> & quNodes,
42 gsMultiPatch<T> & mp)
43 {
44 md.points = quNodes;
45 actives.setZero(7, 1);
46 actives << 0, 1, 2, 3, 4, 5, 6;
47 numActive = actives.rows();
48
49 basisData.setZero(numActive * numActive, md.points.cols());
50 rhsVals.setZero(numActive, md.points.cols());
51
52 gsGeometry<T> & FR = mp.patch(0);
53 gsGeometry<T> & FL = mp.patch(1);
54
55 gsMatrix<T> pointV(FR.parDim(), md.points.cols());
56 pointV.setZero();
57 pointV.row(1) = md.points;
58
59 gsMatrix<T> pointU(FL.parDim(), md.points.cols());
60 pointU.setZero();
61 pointU.row(0) = md.points;
62
63 gsMatrix<T> ones(1, md.points.cols());
64 ones.setOnes();
65
66 gsMatrix<T> lam = ones / math::pow(10,math::ceil(REAL_DIG/2)); // penalty to enforce that alphaL and alphaR will be close to 1. THe penalty depends on the numerical precision
67
68 gsMatrix<T> DuFR(FR.targetDim(), md.points.cols());
69 gsMatrix<T> DvFR(FR.targetDim(), md.points.cols());
70 gsMatrix<T> DvFL(FR.targetDim(), md.points.cols());
71
72 gsMatrix<T> DuFR_DuFR(1, md.points.cols());
73 gsMatrix<T> DvFL_DvFL(1, md.points.cols());
74 gsMatrix<T> DvFR_DvFR(1, md.points.cols());
75
76 gsMatrix<T> DvFL_DuFR(1, md.points.cols());
77 gsMatrix<T> DvFR_DuFR(1, md.points.cols());
78 gsMatrix<T> DvFL_DvFR(1, md.points.cols());
79
80 for(index_t i = 0; i < md.points.cols(); i++)
81 {
82 DuFR.col(i) = FR.jacobian(pointV.col(i)).col(0);
83 DvFR.col(i) = FR.jacobian(pointV.col(i)).col(1); // Same as DuFL
84 DvFL.col(i) = FL.jacobian(pointU.col(i)).col(1);
85
86 // Set scalar product of the jacobian vectors of the geometric mapping
87 DuFR_DuFR.col(i) = DuFR.col(i).transpose() * DuFR.col(i);
88 DvFL_DvFL.col(i) = DvFL.col(i).transpose() * DvFL.col(i);
89 DvFR_DvFR.col(i) = DvFR.col(i).transpose() * DvFR.col(i);
90
91 DvFL_DuFR.col(i) = DvFL.col(i).transpose() * DuFR.col(i);
92 DvFR_DuFR.col(i) = DvFR.col(i).transpose() * DuFR.col(i);
93 DvFL_DvFR.col(i) = DvFL.col(i).transpose() * DvFR.col(i);
94 }
95
96
97
98 // Set Matrix 7x7 values
99// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
100 // alpha_0R
101 basisData.row(0) = 2 * ( md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DvFL - lam) - 2 * md.points.cwiseProduct(DvFL_DvFL - lam) + ones.cwiseProduct(DvFL_DvFL - lam));
102 basisData.row(1) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFL - lam);
103 basisData.row(7) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFL - lam);
104 basisData.row(2) = 2 * (md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DuFR) - 2 * md.points.cwiseProduct(DvFL_DuFR) + ones.cwiseProduct(DvFL_DuFR));
105 basisData.row(14) = 2 * (md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DuFR) - 2 * md.points.cwiseProduct(DvFL_DuFR) + ones.cwiseProduct(DvFL_DuFR));
106 basisData.row(3) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DvFL_DuFR);
107 basisData.row(21) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DvFL_DuFR);
108 basisData.row(4) = 2 * (ones.cwiseProduct(DvFL_DvFR) - 3 * md.points.cwiseProduct(DvFL_DvFR) + 3 * md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DvFR)
109 - md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFL_DvFR));
110 basisData.row(28) = 2 * (ones.cwiseProduct(DvFL_DvFR) - 3 * md.points.cwiseProduct(DvFL_DvFR) + 3 * md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DvFR)
111 - md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFL_DvFR));
112 basisData.row(5) = 4 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
113 basisData.row(35) = 4 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
114 basisData.row(6) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
115 basisData.row(42) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
116
117// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
118 // alpha_1R
119 basisData.row(8) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DvFL - lam);
120 basisData.row(9) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DvFL_DuFR);
121 basisData.row(15) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DvFL_DuFR);
122 basisData.row(10) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DuFR);
123 basisData.row(22) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(DvFL_DuFR);
124 basisData.row(11) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
125 basisData.row(29) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
126 basisData.row(12) = 4 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
127 basisData.row(36) = 4 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFL_DvFR);
128 basisData.row(13) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFL_DvFR);
129 basisData.row(43) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFL_DvFR);
130
131// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
132 // alpha_0L
133 basisData.row(16) = 2 * (md.points.cwiseProduct(md.points).cwiseProduct(DuFR_DuFR - lam) - 2 * md.points.cwiseProduct(DuFR_DuFR - lam) + ones.cwiseProduct(DuFR_DuFR - lam));
134 basisData.row(17) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DuFR_DuFR - lam);
135 basisData.row(23) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(DuFR_DuFR - lam);
136 basisData.row(18) = 2 * (ones.cwiseProduct(DvFR_DuFR) - 3 * md.points.cwiseProduct(DvFR_DuFR) + 3 * md.points.cwiseProduct(md.points).cwiseProduct(DvFR_DuFR)
137 - md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DuFR));
138 basisData.row(30) = 2 * (ones.cwiseProduct(DvFR_DuFR) - 3 * md.points.cwiseProduct(DvFR_DuFR) + 3 * md.points.cwiseProduct(md.points).cwiseProduct(DvFR_DuFR)
139 - md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DuFR));
140 basisData.row(19) = 4 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
141 basisData.row(37) = 4 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
142 basisData.row(20) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
143 basisData.row(44) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
144
145// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
146 // alpha_1L
147 basisData.row(24) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(DuFR_DuFR - lam);
148 basisData.row(25) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
149 basisData.row(31) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
150 basisData.row(26) = 4 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
151 basisData.row(38) = 4 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DuFR);
152 basisData.row(27) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DuFR);
153 basisData.row(45) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DuFR);
154
155// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
156 // beta_0
157 basisData.row(32) = 2 * ( ones.cwiseProduct(DvFR_DvFR) - 4 * md.points.cwiseProduct(DvFR_DvFR) + 6 * md.points.cwiseProduct(md.points).cwiseProduct(DvFR_DvFR)
158 - 4 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DvFR)
159 + md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DvFR));
160 basisData.row(33) = 4 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
161 basisData.row(39) = 4 * md.points.cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
162 basisData.row(34) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
163 basisData.row(46) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
164
165// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
166 // beta_1
167 basisData.row(40) = 8 * md.points.cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
168 basisData.row(41) = 4 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
169 basisData.row(47) = 4 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(ones - md.points).cwiseProduct(DvFR_DvFR);
170
171// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
172 // beta_2
173 basisData.row(48) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(md.points).cwiseProduct(DvFR_DvFR);
174
175// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
176// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
177
178 // Set RHS 7x1 values
179 rhsVals.row(0) = -2 * ( lam - md.points.cwiseProduct(lam) );
180 rhsVals.row(1) = -2 * md.points.cwiseProduct(lam);
181 rhsVals.row(2) = -2 * ( lam - md.points.cwiseProduct(lam) );
182 rhsVals.row(3) = -2 * md.points.cwiseProduct(lam);
183 rhsVals.row(4) = 0 * ones;
184 rhsVals.row(5) = 0 * ones;
185 rhsVals.row(6) = 0 * ones;
186
187 // Initialize local matrix/rhs
188 localMat.setZero(numActive, numActive);
189 localRhs.setZero(numActive, 1 ); //multiple right-hand sides
190 }
191
192
193 // Evaluate on element.
194 inline void evaluateBeta(gsMatrix<T> & quNodes,
195 gsMultiPatch<T> & mp,
196 gsMatrix<T> sol)
197 {
198 GISMO_UNUSED(mp);
199
200 md.points = quNodes;
201 activesBeta.setZero(4, 1);
202 activesBeta << 0, 1, 2, 3;
203 numActiveBeta = activesBeta.rows();
204
205 gsMatrix<T> ones(1, md.points.cols());
206 ones.setOnes();
207
208 gsMatrix<T> alpha_R = sol.row(0) * ( ones - md.points ) + sol.row(1) * md.points;
209 gsMatrix<T> alpha_L = sol.row(2) * ( ones - md.points ) + sol.row(3) * md.points;
210 gsMatrix<T> beta = sol.row(4) * ( md.points.cwiseProduct(md.points) - 2 * md.points + ones ) + 2 * sol.row(5) * md.points.cwiseProduct( ones - md.points ) + sol.row(6) * md.points.cwiseProduct(md.points);
211
212
213 gsMatrix<T> alpha_R_Squared = alpha_R.cwiseProduct(alpha_R);
214 gsMatrix<T> alpha_L_Squared = alpha_L.cwiseProduct(alpha_L);
215
216 gsMatrix<T> alpha_R_L = alpha_R.cwiseProduct(alpha_L);
217
218 gsMatrix<T> lamB = ones / math::pow(10,math::ceil(REAL_DIG/2)); // penalty to enforce that alphaL and alphaR will be close to 1. THe penalty depends on the numerical precision
219
220 basisDataBeta.setZero(numActiveBeta * numActiveBeta, md.points.cols());
221 rhsValsBeta.setZero(numActiveBeta, md.points.cols());
222
223 // Set Matrix 4x4 values
224// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
225 // beta_0R
226 basisDataBeta.row(0) = 2 * ( md.points.cwiseProduct(md.points).cwiseProduct(alpha_L_Squared - lamB) - 2 * md.points.cwiseProduct(alpha_L_Squared - lamB) + ones.cwiseProduct(alpha_L_Squared - lamB));
227 basisDataBeta.row(1) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_L_Squared - lamB);
228 basisDataBeta.row(4) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_L_Squared - lamB);
229 basisDataBeta.row(2) = 2 * ( md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_L) - 2 * md.points.cwiseProduct(alpha_R_L) + ones.cwiseProduct(alpha_R_L));
230 basisDataBeta.row(8) = 2 * ( md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_L) - 2 * md.points.cwiseProduct(alpha_R_L) + ones.cwiseProduct(alpha_R_L));
231 basisDataBeta.row(3) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
232 basisDataBeta.row(12) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
233// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
234 // beta_1R
235 basisDataBeta.row(5) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_L_Squared - lamB);
236 basisDataBeta.row(6) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
237 basisDataBeta.row(9) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
238 basisDataBeta.row(7) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_L);
239 basisDataBeta.row(13) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_L);
240
241
242// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
243 // beta_0L
244 basisDataBeta.row(10) = 2 * ( md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_Squared - lamB) - 2 * md.points.cwiseProduct(alpha_R_Squared - lamB) + ones.cwiseProduct(alpha_R_Squared - lamB));
245 basisDataBeta.row(11) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_Squared - lamB);
246 basisDataBeta.row(14) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_Squared - lamB);
247
248// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
249 // beta_1L
250 basisDataBeta.row(15) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_Squared - lamB);
251
252// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
253// ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
254
255 // Set RHS 4x1 values
256 rhsValsBeta.row(0) = 2 * beta.cwiseProduct(ones - md.points).cwiseProduct(alpha_L);
257 rhsValsBeta.row(1) = 2 * beta.cwiseProduct(md.points).cwiseProduct(alpha_L);
258 rhsValsBeta.row(2) = 2 * beta.cwiseProduct(ones - md.points).cwiseProduct(alpha_R);
259 rhsValsBeta.row(3) = 2 * beta.cwiseProduct(md.points).cwiseProduct(alpha_R);
260
261 // Initialize local matrix/rhs for beta
262 localMatBeta.setZero(numActiveBeta, numActiveBeta);
263 localRhsBeta.setZero(numActiveBeta, 1 ); //multiple right-hand sides
264 }
265
266 inline void assemble(gsDomainIterator<T> & element,
267 const gsVector<T> & quWeights)
268 {
269 GISMO_UNUSED(element);
270
271 for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
272 {
273 // Multiply weight by the geometry measure
274 const T weight = quWeights[k];
275
276 for(index_t i = 0; i < numActive; i++)
277 {
278 localRhs(i, 0) += weight * rhsVals(i,k);
279
280 for(index_t j = 0; j < numActive; j++)
281 {
282 localMat(i, j) += weight * basisData(i*7 + j, k);
283 }
284 }
285 }
286 }
287
288 inline void assembleBeta(gsDomainIterator<T> & element,
289 const gsVector<T> & quWeights)
290 {
291 GISMO_UNUSED(element);
292
293 for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
294 {
295
296 // Multiply weight by the geometry measure
297 const T weight = quWeights[k];
298
299 for(index_t i = 0; i < numActiveBeta; i++)
300 {
301 localRhsBeta(i, 0) += weight * rhsValsBeta(i,k);
302
303 for(index_t j = 0; j < numActiveBeta; j++)
304 {
305 localMatBeta(i, j) += weight * basisDataBeta(i*numActiveBeta + j, k);
306 }
307 }
308 }
309 }
310
311
312 inline void localToGlobal(const gsMatrix<T> & eliminatedDofs,
313 gsSparseSystem<T> & system)
314 {
315 gsMatrix<index_t> actives_temp;
316
317 // Map patch-local DoFs to global DoFs
318 system.mapColIndices(actives, 0, actives_temp);
319 // Add contributions to the system matrix and right-hand side
320 system.push(localMat, localRhs, actives_temp, eliminatedDofs, 0, 0);
321
322 }
323
324
325 inline void localToGlobalBeta(const gsMatrix<T> & eliminatedDofs,
326 gsSparseSystem<T> & system)
327 {
328 gsMatrix<index_t> actives_temp;
329
330 // Map patch-local DoFs to global DoFs
331 system.mapColIndices(activesBeta, 0, actives_temp);
332 // Add contributions to the system matrix and right-hand side
333 system.push(localMatBeta, localRhsBeta, actives_temp, eliminatedDofs, 0, 0);
334
335 }
336
337protected:
338 gsMatrix<index_t> actives;
339 gsMatrix<index_t> activesBeta;
340 gsMatrix<T> basisData;
341 gsMatrix<T> basisDataBeta;
342 index_t numActive;
343 index_t numActiveBeta;
344
345
346protected:
347 // Local values of the right hand side
348 gsMatrix<T> rhsVals;
349 gsMatrix<T> rhsValsBeta;
350
351protected:
352 // Local matrices
353 gsMatrix<T> localMat;
354 gsMatrix<T> localRhs;
355
356 // Local matrices
357 gsMatrix<T> localMatBeta;
358 gsMatrix<T> localRhsBeta;
359
360 gsMapData<T> md;
361
362}; // class gsVisitorGluingData
363
364} // namespace gismo
index_t
#define index_t
Definition gsConfig.h:32
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
gismo
The G+Smo namespace, containing all definitions for the library.