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