G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsC1SurfGluingDataVisitor.h
Go to the documentation of this file.
1 
14 #pragma once
15 
16 namespace gismo
17 {
18 
19 
20 template <class T>
21 class gsC1SurfGluingDataVisitor
22 {
23 public:
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  md.points = quNodes;
199  activesBeta.setZero(4, 1);
200  activesBeta << 0, 1, 2, 3;
201  numActiveBeta = activesBeta.rows();
202 
203  gsMatrix<T> ones(1, md.points.cols());
204  ones.setOnes();
205 
206  gsMatrix<T> alpha_R = sol.row(0) * ( ones - md.points ) + sol.row(1) * md.points;
207  gsMatrix<T> alpha_L = sol.row(2) * ( ones - md.points ) + sol.row(3) * md.points;
208  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);
209 
210 
211  gsMatrix<T> alpha_R_Squared = alpha_R.cwiseProduct(alpha_R);
212  gsMatrix<T> alpha_L_Squared = alpha_L.cwiseProduct(alpha_L);
213 
214  gsMatrix<T> alpha_R_L = alpha_R.cwiseProduct(alpha_L);
215 
216  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
217 
218  basisDataBeta.setZero(numActiveBeta * numActiveBeta, md.points.cols());
219  rhsValsBeta.setZero(numActiveBeta, md.points.cols());
220 
221  // Set Matrix 4x4 values
222 // ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
223  // beta_0R
224  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));
225  basisDataBeta.row(1) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_L_Squared - lamB);
226  basisDataBeta.row(4) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_L_Squared - lamB);
227  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));
228  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));
229  basisDataBeta.row(3) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
230  basisDataBeta.row(12) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
231 // ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
232  // beta_1R
233  basisDataBeta.row(5) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_L_Squared - lamB);
234  basisDataBeta.row(6) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
235  basisDataBeta.row(9) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_L);
236  basisDataBeta.row(7) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_L);
237  basisDataBeta.row(13) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_L);
238 
239 
240 // ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
241  // beta_0L
242  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));
243  basisDataBeta.row(11) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_Squared - lamB);
244  basisDataBeta.row(14) = 2 * md.points.cwiseProduct(ones - md.points).cwiseProduct(alpha_R_Squared - lamB);
245 
246 // ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
247  // beta_1L
248  basisDataBeta.row(15) = 2 * md.points.cwiseProduct(md.points).cwiseProduct(alpha_R_Squared - lamB);
249 
250 // ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
251 // ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
252 
253  // Set RHS 4x1 values
254  rhsValsBeta.row(0) = 2 * beta.cwiseProduct(ones - md.points).cwiseProduct(alpha_L);
255  rhsValsBeta.row(1) = 2 * beta.cwiseProduct(md.points).cwiseProduct(alpha_L);
256  rhsValsBeta.row(2) = 2 * beta.cwiseProduct(ones - md.points).cwiseProduct(alpha_R);
257  rhsValsBeta.row(3) = 2 * beta.cwiseProduct(md.points).cwiseProduct(alpha_R);
258 
259  // Initialize local matrix/rhs for beta
260  localMatBeta.setZero(numActiveBeta, numActiveBeta);
261  localRhsBeta.setZero(numActiveBeta, 1 ); //multiple right-hand sides
262  }
263 
264  inline void assemble(gsDomainIterator<T> & element,
265  const gsVector<T> & quWeights)
266  {
267 
268  for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
269  {
270  // Multiply weight by the geometry measure
271  const T weight = quWeights[k];
272 
273  for(index_t i = 0; i < numActive; i++)
274  {
275  localRhs(i, 0) += weight * rhsVals(i,k);
276 
277  for(index_t j = 0; j < numActive; j++)
278  {
279  localMat(i, j) += weight * basisData(i*7 + j, k);
280  }
281  }
282  }
283  }
284 
285  inline void assembleBeta(gsDomainIterator<T> & element,
286  const gsVector<T> & quWeights)
287  {
288  for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
289  {
290 
291  // Multiply weight by the geometry measure
292  const T weight = quWeights[k];
293 
294  for(index_t i = 0; i < numActiveBeta; i++)
295  {
296  localRhsBeta(i, 0) += weight * rhsValsBeta(i,k);
297 
298  for(index_t j = 0; j < numActiveBeta; j++)
299  {
300  localMatBeta(i, j) += weight * basisDataBeta(i*numActiveBeta + j, k);
301  }
302  }
303  }
304  }
305 
306 
307  inline void localToGlobal(const gsMatrix<T> & eliminatedDofs,
308  gsSparseSystem<T> & system)
309  {
310  gsMatrix<index_t> actives_temp;
311 
312  // Map patch-local DoFs to global DoFs
313  system.mapColIndices(actives, 0, actives_temp);
314  // Add contributions to the system matrix and right-hand side
315  system.push(localMat, localRhs, actives_temp, eliminatedDofs, 0, 0);
316 
317  }
318 
319 
320  inline void localToGlobalBeta(const gsMatrix<T> & eliminatedDofs,
321  gsSparseSystem<T> & system)
322  {
323  gsMatrix<index_t> actives_temp;
324 
325  // Map patch-local DoFs to global DoFs
326  system.mapColIndices(activesBeta, 0, actives_temp);
327  // Add contributions to the system matrix and right-hand side
328  system.push(localMatBeta, localRhsBeta, actives_temp, eliminatedDofs, 0, 0);
329 
330  }
331 
332 protected:
333  gsMatrix<index_t> actives;
334  gsMatrix<index_t> activesBeta;
335  gsMatrix<T> basisData;
336  gsMatrix<T> basisDataBeta;
337  index_t numActive;
338  index_t numActiveBeta;
339 
340 
341 protected:
342  // Local values of the right hand side
343  gsMatrix<T> rhsVals;
344  gsMatrix<T> rhsValsBeta;
345 
346 protected:
347  // Local matrices
348  gsMatrix<T> localMat;
349  gsMatrix<T> localRhs;
350 
351  // Local matrices
352  gsMatrix<T> localMatBeta;
353  gsMatrix<T> localRhsBeta;
354 
355  gsMapData<T> md;
356 
357 }; // class gsVisitorGluingData
358 
359 } // namespace gismo
index_t
#define index_t
Definition: gsConfig.h:32