G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsC1SurfVisitorBasisVertex.h
1
14#pragma once
15
16
17#include <gsUnstructuredSplines/src/gsC1SurfGluingData.h>
18
19
20
21
22namespace gismo
23{
24 template <class T>
25 class gsG1ASVisitorBasisVertex
26 {
27 public:
28
29 gsG1ASVisitorBasisVertex()
30 {
31 }
32
33 void initialize(const gsBasis<T> & basis, //
34 gsQuadRule<T> & rule)
35 {
36 gsVector<index_t> numQuadNodes( basis.dim() );
37 for (index_t i = 0; i < basis.dim(); ++i) // to do: improve
38 numQuadNodes[i] = basis.degree(i) + 1;
39
40 // Setup Quadrature
41 rule = gsGaussRule<T>(numQuadNodes);// NB!
42
43 // Set Geometry evaluation flags
44 md.flags = NEED_MEASURE ;
45
46 localMat.resize(6);
47 localRhs.resize(6);
48
49 rhsVals.resize(6);
50 }
51
52 // Evaluate on element.
53 inline void evaluate(gsBasis<T> & basis, //
54 gsBasis<T> & basis_geo,
55 std::vector<gsBSplineBasis<T>> & basis_plus,
56 std::vector<gsBSplineBasis<T>> & basis_minus,
57 const gsGeometry<T> & geo, // patch
58 gsMatrix<T> & quNodes,
59 gsMatrix<T> & gluingData,
60 std::vector<bool> & isBoundary,
61 gsMatrix<T> & Phi)
62 {
63 GISMO_UNUSED(isBoundary);
64
65 md.points = quNodes;
66
67 // Compute the active basis functions
68 // Assumes actives are the same for all quadrature points on the elements
69 basis.active_into(md.points.col(0), actives);
70
71 // Evaluate basis functions on element
72 basis.eval_into(md.points, basisData);
73
74 // Compute geometry related values
75 geo.computeMap(md);
76
77
78 numActive = actives.rows();
79
80 // Computing c, c+ and c-
81 std::vector<gsMatrix<T>> c_0, c_1;
82 std::vector<gsMatrix < >> c_0_plus, c_1_plus, c_2_plus;
83 std::vector<gsMatrix < >> c_0_plus_deriv, c_1_plus_deriv, c_2_plus_deriv;
84 std::vector<gsMatrix < >> c_0_minus, c_1_minus;
85 for (index_t i = 0; i < 2; i++) // i == 0 == u , i == 1 == v
86 {
87 gsMatrix<T> b_0, b_1;
88 gsMatrix<T> b_0_plus, b_1_plus, b_2_plus;
89 gsMatrix<T> b_0_plus_deriv, b_1_plus_deriv, b_2_plus_deriv;
90 gsMatrix<T> b_0_minus, b_1_minus;
91
92// gsBSplineBasis<T> bsp_temp = dynamic_cast<gsBSplineBasis<T> & >(basis_geo.component(i));
93// T p = bsp_temp.maxDegree();
94// T h_geo = bsp_temp.knots().at(p + 2);
95
96 basis_geo.component(i).evalSingle_into(0, md.points.row(i),b_0); // first
97 basis_geo.component(i).evalSingle_into(1, md.points.row(i),b_1); // second
98
99 basis_plus[i].evalSingle_into(0, md.points.row(i),b_0_plus);
100 basis_plus[i].evalSingle_into(1, md.points.row(i),b_1_plus);
101 basis_plus[i].evalSingle_into(2, md.points.row(i),b_2_plus);
102
103 basis_plus[i].derivSingle_into(0, md.points.row(i),b_0_plus_deriv);
104 basis_plus[i].derivSingle_into(1, md.points.row(i),b_1_plus_deriv);
105 basis_plus[i].derivSingle_into(2, md.points.row(i),b_2_plus_deriv);
106
107 basis_minus[i].evalSingle_into(0, md.points.row(i),b_0_minus);
108 basis_minus[i].evalSingle_into(1, md.points.row(i),b_1_minus);
109
110 // Point zero
111 gsMatrix<T> zero;
112 zero.setZero(2,1);
113
114 gsMatrix<T> b_1_0, b_1_minus_0;
115 basis_geo.component(i).derivSingle_into(1, zero.row(i),b_1_0);
116 basis_minus[i].derivSingle_into(1, zero.row(i),b_1_minus_0);
117// c_0.push_back(b_0 + b_1);
118// c_1.push_back((h_geo / p) * b_1);
119//
120// c_0_minus.push_back(b_0_minus + b_1_minus);
121// c_1_minus.push_back(h_geo/ (p-1) * b_1_minus);
122 T factor_b_1 = 1.0/b_1_0(0,0);
123 c_0.push_back(b_0 + b_1);
124 c_1.push_back(factor_b_1 * b_1);
125
126 T factor_b_1_minus = 1.0/b_1_minus_0(0,0);
127 c_0_minus.push_back(b_0_minus + b_1_minus);
128 c_1_minus.push_back(factor_b_1_minus * b_1_minus);
129
130 gsMatrix<T> der_b_1_plus_0, der2_b_1_plus_0, der2_b_2_plus_0;
131 basis_plus[i].derivSingle_into(1, zero.row(i), der_b_1_plus_0);
132 basis_plus[i].deriv2Single_into(1, zero.row(i), der2_b_1_plus_0);
133 basis_plus[i].deriv2Single_into(2, zero.row(i), der2_b_2_plus_0);
134
135 T factor_c_1_plus = 1/der_b_1_plus_0(0,0);
136 T factor2_c_1_plus = -der2_b_1_plus_0(0,0)/(der_b_1_plus_0(0,0)*der2_b_2_plus_0(0,0));
137 T factor_c_2_plus = 1/der2_b_2_plus_0(0,0);
138
139 c_0_plus.push_back(b_0_plus + b_1_plus + b_2_plus);
140 c_1_plus.push_back(factor_c_1_plus * b_1_plus + factor2_c_1_plus * b_2_plus);
141 c_2_plus.push_back(factor_c_2_plus * b_2_plus );
142
143 c_0_plus_deriv.push_back(b_0_plus_deriv + b_1_plus_deriv + b_2_plus_deriv);
144 c_1_plus_deriv.push_back(factor_c_1_plus * b_1_plus_deriv + factor2_c_1_plus * b_2_plus_deriv);
145 c_2_plus_deriv.push_back(factor_c_2_plus * b_2_plus_deriv);
146 }
147
148// if (g1OptionList.getInt("gluingData") == gluingData::global)
149// {
150// alpha.push_back(gluingData[0].get_alpha_tilde().eval(md.points.row(0))); // u
151// alpha.push_back(gluingData[1].get_alpha_tilde().eval(md.points.row(1))); // v
152// alpha_0.push_back(gluingData[0].get_alpha_tilde().eval(zero.row(0))); // u
153// alpha_0.push_back(gluingData[1].get_alpha_tilde().eval(zero.row(0))); // v
154// alpha_deriv.push_back(gluingData[0].get_alpha_tilde().deriv(zero.row(0))); // u
155// alpha_deriv.push_back(gluingData[1].get_alpha_tilde().deriv(zero.row(0))); // v
156//
157// beta.push_back(gluingData[0].get_beta_tilde().eval(md.points.row(0))); // u
158// beta.push_back(gluingData[1].get_beta_tilde().eval(md.points.row(1))); // v
159// beta_0.push_back(gluingData[0].get_beta_tilde().eval(zero.row(0))); // u
160// beta_0.push_back(gluingData[1].get_beta_tilde().eval(zero.row(0))); // v
161// beta_deriv.push_back(gluingData[0].get_beta_tilde().deriv(zero.row(0))); // u
162// beta_deriv.push_back(gluingData[1].get_beta_tilde().deriv(zero.row(0))); // v
163//
164// }
165// else if (g1OptionList.getInt("gluingData") == gluingData::local)
166// {
167// alpha.push_back(gluingData[0].get_local_alpha_tilde(0).eval(md.points.row(0))); // u
168// alpha.push_back(gluingData[1].get_local_alpha_tilde(0).eval(md.points.row(1))); // v
169// alpha_0.push_back(gluingData[0].get_local_alpha_tilde(0).eval(zero.row(0))); // u
170// alpha_0.push_back(gluingData[1].get_local_alpha_tilde(0).eval(zero.row(0))); // v
171// alpha_deriv.push_back(gluingData[0].get_local_alpha_tilde(0).deriv(zero.row(0))); // u
172// alpha_deriv.push_back(gluingData[1].get_local_alpha_tilde(0).deriv(zero.row(0))); // v
173//
174// beta.push_back(gluingData[0].get_local_beta_tilde(0).eval(md.points.row(0))); // u
175// beta.push_back(gluingData[1].get_local_beta_tilde(0).eval(md.points.row(1))); // v
176// beta_0.push_back(gluingData[0].get_local_beta_tilde(0).eval(zero.row(0))); // u
177// beta_0.push_back(gluingData[1].get_local_beta_tilde(0).eval(zero.row(0))); // v
178// beta_deriv.push_back(gluingData[0].get_local_beta_tilde(0).deriv(zero.row(0))); // u
179// beta_deriv.push_back(gluingData[1].get_local_beta_tilde(0).deriv(zero.row(0))); // v
180// }
181
182 std::vector<gsMatrix<T>> alpha, beta, alpha_0, beta_0, alpha_deriv, beta_deriv;
183
184 // Point zero
185 gsMatrix<T> zero;
186 zero.setZero(2,1);
187
188 // Geo data:
189 gsMatrix<T> geo_jac, geo_der2;
190 geo_jac = geo.jacobian(zero);
191 geo_der2 = geo.deriv2(zero);
192
193 gsMatrix<T> zeros(1, md.points.cols());
194 zeros.setZero();
195
196 // Point One
197 gsMatrix<T> one;
198 one.setOnes(2,1);
199
200 gsMatrix<T> ones(1, md.points.cols());
201 ones.setOnes();
202
203 alpha.push_back( gluingData(0, 0) * ( ones - md.points.row(0) ) + gluingData(1, 0) * md.points.row(0) ); // u
204 alpha.push_back( gluingData(0, 1) * ( ones - md.points.row(1) ) + gluingData(1, 1) * md.points.row(1) ); // v
205
206
207 alpha_0.push_back( gluingData(0, 0) * ( one.row(0) - zero.row(0) ) + gluingData(1, 0) * zero.row(0) ); // u
208 alpha_0.push_back( gluingData(0, 1) * ( one.row(0) - zero.row(0) ) + gluingData(1, 1) * zero.row(0) ); // v
209 alpha_deriv.push_back( ( gluingData(1, 0) - gluingData(0, 0) ) * ones.col(0) ); // u
210 alpha_deriv.push_back( ( gluingData(1, 1) - gluingData(0, 1) ) * ones.col(0) ); // v
211
212 beta.push_back( gluingData(2, 0) * ( ones - md.points.row(0) ) + gluingData(3, 0) * md.points.row(0) ); // u
213 beta.push_back( gluingData(2, 1) * ( ones - md.points.row(1) ) + gluingData(3, 1) * md.points.row(1) ); // v
214 beta_0.push_back( gluingData(2, 0) * ( one.row(0) - zero.row(0) ) + gluingData(3, 0) * zero.row(0) ); // u
215 beta_0.push_back( gluingData(2, 1) * ( one.row(0) - zero.row(0) ) + gluingData(3, 1) * zero.row(0) ); // v
216 beta_deriv.push_back( ( gluingData(3, 0) - gluingData(2, 0) ) * ones.col(0) ); // u
217 beta_deriv.push_back( ( gluingData(3, 1) - gluingData(2, 1) ) * ones.col(0) ); // v
218
219 // Compute dd^^(i_k) and dd^^(i_k-1)
220 gsMatrix<T> dd_ik_plus, dd_ik_minus;
221 gsMatrix<T> dd_ik_minus_deriv, dd_ik_plus_deriv;
222
223 dd_ik_minus = ( -1 / alpha_0[0](0,0) ) * ( geo_jac.col(1) +
224 beta_0[0](0,0) * geo_jac.col(0) );
225
226 dd_ik_plus = ( 1 / alpha_0[1](0,0) ) * ( geo_jac.col(0) +
227 beta_0[1](0,0) * geo_jac.col(1) );
228
229 gsMatrix<T> geo_deriv2_12(geo.targetDim(), 1), geo_deriv2_11(geo.targetDim(), 1), geo_deriv2_22(geo.targetDim(), 1);
230
231 geo_deriv2_12.row(0) = geo_der2.row(2);
232 geo_deriv2_12.row(1) = geo_der2.row(5);
233
234 geo_deriv2_11.row(0) = geo_der2.row(0);
235 geo_deriv2_11.row(1) = geo_der2.row(3);
236
237 geo_deriv2_22.row(0) = geo_der2.row(1);
238 geo_deriv2_22.row(1) = geo_der2.row(4);
239
240 if(geo.parDim() + 1 == geo.targetDim()) // In the surface case the dimension of the second derivative vector is 9x1
241 {
242 geo_deriv2_12.row(2) = geo_der2.row(8);
243
244 geo_deriv2_11.row(2) = geo_der2.row(6);
245
246 geo_deriv2_22.row(2) = geo_der2.row(7);
247 }
248
249 gsMatrix<T> alpha_squared_u = alpha_0[0]*alpha_0[0];
250 gsMatrix<T> alpha_squared_v = alpha_0[1]*alpha_0[1];
251
252 dd_ik_minus_deriv = -1/(alpha_squared_u(0,0)) * // N^2
253 ( ( geo_deriv2_12 + (beta_deriv[0](0,0) * geo_jac.col(0) + beta_0[0](0,0) * geo_deriv2_11) ) * alpha_0[0](0,0) -
254 ( geo_jac.col(1) + beta_0[0](0,0) * geo_jac.col(0) ) * alpha_deriv[0](0,0) );
255
256
257 dd_ik_plus_deriv = 1/(alpha_squared_v(0,0)) *
258 ( ( geo_deriv2_12 + (beta_deriv[1](0,0) * geo_jac.col(1) + beta_0[1](0,0) * geo_deriv2_22) ) * alpha_0[1](0,0) -
259 ( geo_jac.col(0) + beta_0[1](0,0) * geo_jac.col(1) ) * alpha_deriv[1](0,0) );
260
261 //if (isBoundary[0] == false)
262 // gsInfo << dd_ik_minus_deriv << "\n";
263 //if (isBoundary[1] == false)
264 // gsInfo << dd_ik_plus_deriv << "\n";
265
266 // Comupute d_(0,0)^(i_k), d_(1,0)^(i_k), d_(0,1)^(i_k), d_(1,1)^(i_k) ; i_k == 2
267 std::vector<gsMatrix<T>> d_ik;
268 d_ik.push_back(Phi.row(0).transpose());
269
270 d_ik.push_back(Phi.block(1, 0, geo.targetDim(), 6).transpose() * geo_jac.col(0) ); // deriv into u
271
272
273// gsInfo << "d_ik: " << d_ik.back() << "\n";
274// gsInfo << "======================================= \n";
275// gsInfo << "geo_jac.col(0): " << geo_jac.col(0) << "\n";
276// gsInfo << "======================================= \n";
277
278
279 d_ik.push_back(Phi.block(1, 0, geo.targetDim(), 6).transpose() * geo_jac.col(1) ); // deriv into v
280
281 // Hessian
282 if(geo.parDim() + 1 == geo.targetDim()) // In the surface case the dimension of the second derivative vector is 9x1
283 {
284 d_ik.push_back( (geo_jac(0,0) * Phi.row(4).transpose() + geo_jac(1,0) * Phi.row(7).transpose() + geo_jac(2,0) * Phi.row(10).transpose()) * geo_jac(0,1) +
285 (geo_jac(0,0) * Phi.row(5).transpose() + geo_jac(1,0) * Phi.row(8).transpose() + geo_jac(2,0) * Phi.row(11).transpose()) * geo_jac(1,1) +
286 (geo_jac(0,0) * Phi.row(6).transpose() + geo_jac(1,0) * Phi.row(9).transpose() + geo_jac(2,0) * Phi.row(12).transpose()) * geo_jac(2,1) +
287 Phi.block(1, 0, 1, 6).transpose() * geo_der2.row(2) +
288 Phi.block(2, 0, 1, 6).transpose() * geo_der2.row(5) +
289 Phi.block(3, 0, 1, 6).transpose() * geo_der2.row(8) );
290
291 }
292 else
293 {
294 d_ik.push_back( (geo_jac(0,0) * Phi.col(3) + geo_jac(1,0) * Phi.col(4)) * geo_jac(0,1) +
295 (geo_jac(0,0) * Phi.col(4) + geo_jac(1,0) * Phi.col(5)) * geo_jac(1,1) +
296 Phi.block(0,1, 6,1) * geo_der2.row(2) +
297 Phi.block(0,2, 6,1) * geo_der2.row(5)); // Hessian
298 }
299
300// gsInfo << "d_ik: " << d_ik.back() << " : " << Phi << "\n";
301 // Compute d_(*,*)^(il,ik)
302 std::vector<gsMatrix<T>> d_ilik_minus, d_ilik_plus;
303
304// d_(*,*)^(ik-1,ik)
305 d_ilik_minus.push_back(Phi.row(0).transpose());
306
307 d_ilik_minus.push_back(Phi.block(1, 0, geo.targetDim(), 6).transpose() * geo_jac.col(0));
308
309 if(geo.parDim() + 1 == geo.targetDim()) // In the surface case the dimension of the second derivative vector is 9x1
310 {
311 d_ilik_minus.push_back( (geo_jac(0,0) * Phi.row(4).transpose() + geo_jac(1,0) * Phi.row(7).transpose() + geo_jac(2,0) * Phi.row(10).transpose()) * geo_jac(0,0) +
312 (geo_jac(0,0) * Phi.row(5).transpose() + geo_jac(1,0) * Phi.row(8).transpose() + geo_jac(2,0) * Phi.row(11).transpose()) * geo_jac(1,0) +
313 (geo_jac(0,0) * Phi.row(6).transpose() + geo_jac(1,0) * Phi.row(9).transpose() + geo_jac(2,0) * Phi.row(12).transpose()) * geo_jac(2,0) +
314 Phi.block(1, 0, 1, 6).transpose() * geo_der2.row(0) +
315 Phi.block(2, 0, 1, 6).transpose() * geo_der2.row(3) +
316 Phi.block(3, 0, 1, 6).transpose() * geo_der2.row(6) );
317 }
318 else
319 {
320 d_ilik_minus.push_back( (geo_jac(0,0) * Phi.col(3) + geo_jac(1,0) * Phi.col(4))*geo_jac(0,0) +
321 (geo_jac(0,0) * Phi.col(4) + geo_jac(1,0) * Phi.col(5))*geo_jac(1,0) +
322 Phi.block(0,1,6,1) * geo_der2.row(0) +
323 Phi.block(0,2,6,1) * geo_der2.row(3));
324 }
325
326 d_ilik_minus.push_back(Phi.block(1, 0, geo.targetDim(), 6).transpose() * dd_ik_minus);
327
328 if(geo.parDim() + 1 == geo.targetDim()) // In the surface case the dimension of the second derivative vector is 9x1
329 {
330 d_ilik_minus.push_back( (geo_jac(0,0) * Phi.row(4).transpose() + geo_jac(1,0) * Phi.row(7).transpose() + geo_jac(2,0) * Phi.row(10).transpose()) * dd_ik_minus(0,0) +
331 (geo_jac(0,0) * Phi.row(5).transpose() + geo_jac(1,0) * Phi.row(8).transpose() + geo_jac(2,0) * Phi.row(11).transpose()) * dd_ik_minus(1,0) +
332 (geo_jac(0,0) * Phi.row(6).transpose() + geo_jac(1,0) * Phi.row(9).transpose() + geo_jac(2,0) * Phi.row(12).transpose()) * dd_ik_minus(2,0) +
333 Phi.block(1, 0, 1, 6).transpose() * dd_ik_minus_deriv.row(0) +
334 Phi.block(2, 0, 1, 6).transpose() * dd_ik_minus_deriv.row(1) +
335 Phi.block(3, 0, 1, 6).transpose() * dd_ik_minus_deriv.row(2) );
336 }
337 else
338 {
339 d_ilik_minus.push_back( (geo_jac(0,0) * Phi.col(3) + geo_jac(1,0) * Phi.col(4)) * dd_ik_minus(0,0) +
340 (geo_jac(0,0) * Phi.col(4) + geo_jac(1,0) * Phi.col(5)) * dd_ik_minus(1,0) +
341 Phi.block(0,1,6,1) * dd_ik_minus_deriv.row(0) +
342 Phi.block(0,2,6,1) * dd_ik_minus_deriv.row(1));
343 }
344
345
346// d_(*,*)^(ik+1,ik)
347 d_ilik_plus.push_back(Phi.row(0).transpose());
348
349 d_ilik_plus.push_back(Phi.block(1, 0, geo.targetDim(), 6).transpose() * geo_jac.col(1));
350
351 if(geo.parDim() + 1 == geo.targetDim()) // In the surface case the dimension of the second derivative vector is 9x1
352 {
353 d_ilik_plus.push_back( (geo_jac(0,1) * Phi.row(4).transpose() + geo_jac(1,1) * Phi.row(7).transpose() + geo_jac(2,1) * Phi.row(10).transpose()) * geo_jac(0,1) +
354 (geo_jac(0,1) * Phi.row(5).transpose() + geo_jac(1,1) * Phi.row(8).transpose() + geo_jac(2,1) * Phi.row(11).transpose()) * geo_jac(1,1) +
355 (geo_jac(0,1) * Phi.row(6).transpose() + geo_jac(1,1) * Phi.row(9).transpose() + geo_jac(2,1) * Phi.row(12).transpose()) * geo_jac(2,1) +
356 Phi.block(1, 0, 1, 6).transpose() * geo_der2.row(1) +
357 Phi.block(2, 0, 1, 6).transpose() * geo_der2.row(4) +
358 Phi.block(3, 0, 1, 6).transpose() * geo_der2.row(7) );
359 }
360 else
361 {
362 d_ilik_plus.push_back( (geo_jac(0,1) * Phi.col(3) + geo_jac(1,1) * Phi.col(4)) * geo_jac(0,1) +
363 (geo_jac(0,1) * Phi.col(4) + geo_jac(1,1) * Phi.col(5)) * geo_jac(1,1) +
364 Phi.block(0,1,6,1) * geo_der2.row(1) +
365 Phi.block(0,2,6,1) * geo_der2.row(4) );
366 }
367
368 d_ilik_plus.push_back(Phi.block(1, 0, geo.targetDim(), 6).transpose() * dd_ik_plus);
369
370 if(geo.parDim() + 1 == geo.targetDim()) // In the surface case the dimension of the second derivative vector is 9x1
371 {
372 d_ilik_plus.push_back( (geo_jac(0,1) * Phi.row(4).transpose() + geo_jac(1,1) * Phi.row(7).transpose() + geo_jac(2,1) * Phi.row(10).transpose()) * dd_ik_plus(0,0) +
373 (geo_jac(0,1) * Phi.row(5).transpose() + geo_jac(1,1) * Phi.row(8).transpose() + geo_jac(2,1) * Phi.row(11).transpose()) * dd_ik_plus(1,0) +
374 (geo_jac(0,1) * Phi.row(6).transpose() + geo_jac(1,1) * Phi.row(9).transpose() + geo_jac(2,1) * Phi.row(12).transpose()) * dd_ik_plus(2,0) +
375 Phi.block(1, 0, 1, 6).transpose() * dd_ik_plus_deriv.row(0) +
376 Phi.block(2, 0, 1, 6).transpose() * dd_ik_plus_deriv.row(1) +
377 Phi.block(3, 0, 1, 6).transpose() * dd_ik_plus_deriv.row(2) );
378 }
379 else
380 {
381 d_ilik_plus.push_back( (geo_jac(0,1) * Phi.col(3) + geo_jac(1,1) * Phi.col(4)) * dd_ik_plus(0,0) +
382 (geo_jac(0,1) * Phi.col(4) + geo_jac(1,1) * Phi.col(5)) * dd_ik_plus(1,0) +
383 Phi.block(0,1,6,1) * dd_ik_plus_deriv.row(0) +
384 Phi.block(0,2,6,1) * dd_ik_plus_deriv.row(1) );
385 }
386
387
388 for (index_t i = 0; i < 6; i++)
389 {
390 rhsVals.at(i) = d_ilik_minus.at(0)(i,0) * (c_0_plus.at(0).cwiseProduct(c_0.at(1)) -
391 beta[0].cwiseProduct(c_0_plus_deriv.at(0).cwiseProduct(c_1.at(1)))) +
392 d_ilik_minus.at(1)(i,0) * (c_1_plus.at(0).cwiseProduct(c_0.at(1)) -
393 beta[0].cwiseProduct(c_1_plus_deriv.at(0).cwiseProduct(c_1.at(1)))) +
394 d_ilik_minus.at(2)(i,0) * (c_2_plus.at(0).cwiseProduct(c_0.at(1)) -
395 beta[0].cwiseProduct(c_2_plus_deriv.at(0).cwiseProduct(c_1.at(1)))) -
396 d_ilik_minus.at(3)(i,0) * alpha[0].cwiseProduct(c_0_minus.at(0).cwiseProduct(c_1.at(1))) -
397 d_ilik_minus.at(4)(i,0) * alpha[0].cwiseProduct(c_1_minus.at(0).cwiseProduct(c_1.at(1))); // f*_(ik-1,ik)
398
399 rhsVals.at(i) += d_ilik_plus.at(0)(i,0) * (c_0_plus.at(1).cwiseProduct(c_0.at(0)) -
400 beta[1].cwiseProduct(c_0_plus_deriv.at(1).cwiseProduct(c_1.at(0)))) +
401 d_ilik_plus.at(1)(i,0) * (c_1_plus.at(1).cwiseProduct(c_0.at(0)) -
402 beta[1].cwiseProduct(c_1_plus_deriv.at(1).cwiseProduct(c_1.at(0)))) +
403 d_ilik_plus.at(2)(i,0) * (c_2_plus.at(1).cwiseProduct(c_0.at(0)) -
404 beta[1].cwiseProduct(c_2_plus_deriv.at(1).cwiseProduct(c_1.at(0)))) +
405 d_ilik_plus.at(3)(i,0) * alpha[1].cwiseProduct(c_0_minus.at(1).cwiseProduct(c_1.at(0))) +
406 d_ilik_plus.at(4)(i,0) * alpha[1].cwiseProduct(c_1_minus.at(1).cwiseProduct(c_1.at(0))); // f*_(ik+1,ik)
407
408 rhsVals.at(i) -= d_ik.at(0)(i,0) * c_0.at(0).cwiseProduct(c_0.at(1)) + d_ik.at(2)(i,0) * c_0.at(0).cwiseProduct(c_1.at(1)) +
409 d_ik.at(1)(i,0) * c_1.at(0).cwiseProduct(c_0.at(1)) + d_ik.at(3)(i,0) * c_1.at(0).cwiseProduct(c_1.at(1)); // f*_(ik)
410
411 localMat.at(i).setZero(numActive, numActive);
412 localRhs.at(i).setZero(numActive, rhsVals.at(i).rows());//multiple right-hand sides
413
414 localMat.at(i).setZero(numActive, numActive);
415 localRhs.at(i).setZero(numActive, rhsVals.at(i).rows());//multiple right-hand sides
416 }
417
418
419
420 } // evaluate
421
422 inline void assemble(gsDomainIterator<T> & element,
423 const gsVector<T> & quWeights)
424 {
425 GISMO_UNUSED(element);
426
427 gsMatrix<T> & basisVals = basisData;
428 for (index_t i = 0; i < 6; i++)
429 {
430 // ( u, v)
431 localMat.at(i).noalias() = basisData * quWeights.asDiagonal() * md.measures.asDiagonal() * basisData.transpose();
432
433 for (index_t k = 0; k < quWeights.rows(); ++k) // loop over quadrature nodes
434 {
435 T weight = quWeights[k];
436 gsMatrix<T> Jk = md.jacobian(k);
437
438 if( Jk.dim().second + 1 == Jk.dim().first)
439 {
440 gsMatrix<T> G = Jk.transpose() * Jk;
441 T detG = G.determinant();
442 weight *= sqrt(detG);
443 }
444 else
445 {
446 weight *= md.measure(k);
447 }
448 // Multiply weight by the geometry measure
449
450 localRhs.at(i).noalias() += weight * (basisVals.col(k) * rhsVals.at(i).col(k).transpose());
451 }
452 }
453 }
454
455 inline void localToGlobal(const index_t patchIndex,
456 const std::vector<gsMatrix<T> > & eliminatedDofs,
457 std::vector< gsSparseSystem<T> > & system)
458 {
459 gsMatrix<index_t> actives_temp;
460 for (size_t i = 0; i < system.size(); i++) // 6
461 {
462 // Map patch-local DoFs to global DoFs
463 system.at(i).mapColIndices(actives, patchIndex, actives_temp);
464 // Add contributions to the system matrix and right-hand side
465 system.at(i).push(localMat.at(i), localRhs.at(i), actives_temp, eliminatedDofs[0], 0, 0);
466 }
467 }
468
469 protected:
470 gsMatrix<index_t> actives;
471 gsMatrix<T> basisData;
472 index_t numActive;
473
474 protected:
475 // Local values of the right hand side
476 std::vector< gsMatrix<T> > rhsVals;
477
478 protected:
479 // Local matrices
480 std::vector< gsMatrix<T> > localMat;
481 std::vector< gsMatrix<T> > localRhs;
482
483 gsMapData<T> md;
484
485 }; // class gsVisitorG1BasisVertex
486} // namespace gismo
gismo::gsBSplineBasis
A univariate B-spline basis.
Definition gsBSplineBasis.h:700
gismo::gsBasis
A basis represents a family of scalar basis functions defined over a common parameter domain.
Definition gsBasis.h:79
gismo::gsBasis::component
virtual const gsBasis< T > & component(short_t i) const
For a tensor product basis, return the (const) 1-d basis for the i-th parameter component.
Definition gsBasis.hpp:563
gismo::gsDomainIterator
Class which enables iteration over all elements of a parameter domain.
Definition gsDomainIterator.h:68
gismo::gsFunctionSet::deriv2
gsMatrix< T > deriv2(const gsMatrix< T > &u) const
Evaluates the second derivatives of active (i.e., non-zero) functions at points u.
Definition gsFunctionSet.hpp:138
gismo::gsFunction::computeMap
virtual void computeMap(gsMapData< T > &InOut) const
Computes map function data.
Definition gsFunction.hpp:817
gismo::gsGaussRule
Class that represents the (tensor) Gauss-Legendre quadrature rule.
Definition gsGaussRule.h:28
gismo::gsGeometry
Abstract base class representing a geometry map.
Definition gsGeometry.h:93
gismo::gsGeometry::targetDim
short_t targetDim() const
Dimension of the ambient physical space (overriding gsFunction::targetDim())
Definition gsGeometry.h:286
gismo::gsGeometry::parDim
short_t parDim() const
Dimension d of the parameter domain (same as domainDim()).
Definition gsGeometry.hpp:190
gismo::gsMapData
the gsMapData is a cache of pre-computed function (map) values.
Definition gsFuncData.h:349
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition gsMatrix.h:211
gismo::gsQuadRule
Class representing a reference quadrature rule.
Definition gsQuadRule.h:29
gismo::gsSparseSystem
A sparse linear system indexed by sets of degrees of freedom.
Definition gsSparseSystem.h:30
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
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.
gismo::NEED_MEASURE
@ NEED_MEASURE
The density of the measure pull back.
Definition gsForwardDeclarations.h:76