Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__CONTINUITY_CONSTRAINTS__P1_CONTINUITY_CONSTRAINTS_IMPL__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__CONTINUITY_CONSTRAINTS__P1_CONTINUITY_CONSTRAINTS_IMPL__
35 :
36 : #include "lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints.h"
37 : #include "lib_disc/domain.h"
38 :
39 : namespace ug {
40 :
41 : /// sets a matrix row corresponding to averaging the constrained index
42 : template <typename TMatrix>
43 0 : void SetInterpolation(TMatrix& A,
44 : std::vector<size_t> & constrainedIndex,
45 : std::vector<std::vector<size_t> >& vConstrainingIndex,
46 : bool assembleLinearProblem = true)
47 : {
48 : //typedef typename TMatrix::row_iterator row_iterator;
49 : //typedef typename TMatrix::value_type block_type;
50 :
51 : // check number of indices passed
52 : for(size_t i = 0; i < vConstrainingIndex.size(); ++i)
53 : UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(),
54 : "Wrong number of indices.");
55 :
56 : // loop all constrained dofs
57 0 : for(size_t i = 0; i < constrainedIndex.size(); ++i)
58 : {
59 0 : const size_t ai = constrainedIndex[i];
60 :
61 : // remove all couplings
62 0 : SetRow(A, ai, 0.0);
63 :
64 : // set diag of row to identity
65 0 : A(ai, ai) = 1.0;
66 :
67 : // set coupling to all constraining dofs the inverse of the
68 : // number of constraining dofs
69 : // This is only required if assembling for a linear problem.
70 0 : if (assembleLinearProblem)
71 : {
72 0 : const number frac = -1.0/(vConstrainingIndex.size());
73 0 : for(size_t j=0; j < vConstrainingIndex.size(); ++j)
74 0 : A(ai, vConstrainingIndex[j][i]) = frac;
75 : }
76 : }
77 0 : }
78 :
79 : template <typename TVector>
80 0 : void InterpolateValues(TVector& u,
81 : std::vector<size_t>& constrainedIndex,
82 : std::vector<std::vector<size_t> >& vConstrainingIndex)
83 : {
84 : typedef typename TVector::value_type block_type;
85 :
86 : // check number of indices passed
87 0 : for(size_t i = 0; i < vConstrainingIndex.size(); ++i)
88 : UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(),
89 : "Wrong number of indices.");
90 :
91 : // scaling factor
92 0 : const number frac = 1./(vConstrainingIndex.size());
93 :
94 : // loop constrained indices
95 0 : for(size_t i = 0; i < constrainedIndex.size(); ++i)
96 : {
97 : // get value to be interpolated
98 0 : block_type& val = u[constrainedIndex[i]];
99 :
100 : // reset value
101 0 : val = 0.0;
102 :
103 : // add equally from all constraining indices
104 0 : for(size_t j=0; j < vConstrainingIndex.size(); ++j)
105 0 : VecScaleAdd(val, 1.0, val, frac, u[vConstrainingIndex[j][i]]);
106 : //val += frac * u[vConstrainingIndex[j][i]];
107 : }
108 0 : }
109 :
110 :
111 :
112 : template <typename TMatrix>
113 0 : void SplitAddRow_Symmetric(TMatrix& A,
114 : std::vector<size_t>& constrainedIndex,
115 : std::vector<std::vector<size_t> >& vConstrainingIndex)
116 : {
117 : typedef typename TMatrix::value_type block_type;
118 : typedef typename TMatrix::row_iterator row_iterator;
119 :
120 : size_t nConstrg = vConstrainingIndex.size();
121 : UG_ASSERT(nConstrg, "There have to be constraining indices!");
122 :
123 : // check number of indices passed
124 : for(size_t i = 0; i < nConstrg; ++i)
125 : UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(),
126 : "Wrong number of indices.");
127 :
128 : // scaling factor
129 0 : const number frac = 1.0 / nConstrg;
130 :
131 : // handle each constrained index
132 0 : for(size_t i = 0; i < constrainedIndex.size(); ++i)
133 : {
134 0 : const size_t algInd = constrainedIndex[i];
135 :
136 : // add coupling constrained index -> constrained index
137 : // we don't have to adjust the block itself, since the row of
138 : // constraints will be set to interpolation afterwards
139 0 : block_type diagEntry = A(algInd, algInd);
140 :
141 : // scale by weight
142 0 : diagEntry *= frac*frac;
143 :
144 : // add coupling constrained dof -> constrained dof
145 0 : for(size_t k = 0; k < vConstrainingIndex.size(); ++k)
146 0 : for(size_t m = 0; m < vConstrainingIndex.size(); ++m)
147 0 : A(vConstrainingIndex[k][i], vConstrainingIndex[m][i]) += diagEntry;
148 :
149 : #if 0
150 :
151 : // (handled separately as it would otherwise trigger an assert -
152 : // manipulation of matrix row while iterating over it)
153 : const block_type block = A(algInd, algInd);
154 : size_t nBlockCols = GetCols(block);
155 : UG_ASSERT(nBlockCols == constrainedIndex.size(),
156 : "Number of block cols and number of constrained DoFs in hanging vertex not equal.");
157 : for (size_t blockIndConn = 0; blockIndConn < nBlockCols; ++blockIndConn)
158 : {
159 : const number val = BlockRef(block, blockInd, blockIndConn) * frac*frac;
160 : for (size_t k = 0; k < nConstrg; ++k)
161 : for (size_t m = 0; m < nConstrg; ++m)
162 : DoFRef(A, vConstrainingIndex[k][i], vConstrainingIndex[m][blockIndConn]) += val;
163 : }
164 : #endif
165 :
166 : // loop coupling between constrained dof -> any other dof
167 0 : for(row_iterator conn = A.begin_row(algInd); conn != A.end_row(algInd); ++conn)
168 : {
169 : const size_t algIndConn = conn.index();
170 :
171 : // diagonal entry already handled
172 0 : if (algIndConn == algInd)
173 0 : continue;
174 :
175 : // warning: do NOT use references here!
176 : // they might become invalid when A is accessed at an entry that does not exist yet
177 : // FIXME: This will only work properly if there is an entry A(i,j) for any entry
178 : // A(j,i) in the matrix! Is this always the case!?
179 0 : block_type block = conn.value();
180 0 : block_type blockT = A(algIndConn, algInd);
181 :
182 : // multiply the cpl value by the inverse number of constraining
183 0 : block *= frac;
184 0 : blockT *= frac;
185 :
186 : // add the coupling to the constraining indices rows
187 0 : for (size_t k = 0; k < nConstrg; ++k)
188 : {
189 : UG_ASSERT(vConstrainingIndex[k][i] != constrainedIndex[i],
190 : "Modifying 'this' (=conn referenced) matrix row is invalid!" << constrainedIndex[i]);
191 0 : A(vConstrainingIndex[k][i], algIndConn) += block;
192 0 : A(algIndConn, vConstrainingIndex[k][i]) += blockT;
193 : }
194 :
195 : // set the split coupling to zero
196 : // this needs to be done only in columns, since the row associated to
197 : // the constrained index will be set to unity (or interpolation).
198 0 : A(algIndConn, algInd) = 0.0;
199 : }
200 : }
201 0 : }
202 :
203 : template <typename TMatrix>
204 0 : void SplitAddRow_OneSide(TMatrix& A,
205 : std::vector<size_t>& constrainedIndex,
206 : std::vector<std::vector<size_t> >& vConstrainingIndex)
207 : {
208 : typedef typename TMatrix::value_type block_type;
209 : typedef typename TMatrix::row_iterator row_iterator;
210 :
211 : UG_ASSERT(!vConstrainingIndex.empty(), "There have to be constraining indices!");
212 :
213 : // check number of indices passed
214 0 : for(size_t i = 0; i < vConstrainingIndex.size(); ++i)
215 : UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(),
216 : "Wrong number of indices.");
217 :
218 : // scaling factor
219 0 : const number frac = 1./(vConstrainingIndex.size());
220 :
221 0 : for(size_t i = 0; i < constrainedIndex.size(); ++i)
222 : {
223 0 : const size_t algInd = constrainedIndex[i];
224 :
225 : // choose randomly the first dof to add whole row
226 0 : const size_t addTo = vConstrainingIndex[0][i];
227 :
228 0 : block_type diagEntry = A(algInd, algInd);
229 :
230 : // scale by weight
231 0 : diagEntry *= frac;
232 :
233 : // add coupling
234 0 : for(size_t k = 0; k < vConstrainingIndex.size(); ++k)
235 0 : A(addTo, vConstrainingIndex[k][i]) += diagEntry;
236 :
237 : // loop coupling between constrained dof -> any other dof
238 0 : for(row_iterator conn = A.begin_row(algInd); conn != A.end_row(algInd); ++conn)
239 : {
240 : const size_t algIndConn = conn.index();
241 :
242 : // skip self-coupling (already handled)
243 0 : if (algIndConn == algInd) continue;
244 :
245 : // warning: do NOT use references here!
246 : // they might become invalid when A is accessed at an entry that does not exist yet
247 : // FIXME: This will only work properly if there is an entry A(i,j) for any entry
248 : // A(j,i) in the matrix! Is this always the case!?
249 0 : const block_type block = conn.value();
250 0 : block_type blockT = A(algIndConn, algInd);
251 :
252 0 : blockT *= frac;
253 :
254 : // add the coupling to the constraining indices rows
255 0 : for(size_t k = 0; k < vConstrainingIndex.size(); ++k)
256 0 : A(algIndConn, vConstrainingIndex[k][i]) += blockT;
257 :
258 : // coupling due to one side adding
259 0 : A(addTo, algIndConn) += block;
260 :
261 : // set the split coupling to zero
262 : // this must only be done in columns, since the row associated to
263 : // the constrained index will be set to an interpolation.
264 0 : A(algIndConn, algInd) = 0.0;
265 : }
266 : }
267 0 : }
268 :
269 : template <typename TVector>
270 0 : void SplitAddRhs_Symmetric(TVector& rhs,
271 : std::vector<size_t> & constrainedIndex,
272 : std::vector<std::vector<size_t> >& vConstrainingIndex)
273 : {
274 : typedef typename TVector::value_type block_type;
275 :
276 : // check number of indices passed
277 0 : for(size_t i = 0; i < vConstrainingIndex.size(); ++i)
278 : UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(),
279 : "Wrong number of indices.");
280 :
281 : // scaling factor
282 0 : const number frac = 1./(vConstrainingIndex.size());
283 :
284 : // loop constrained indices
285 0 : for(size_t i = 0; i < constrainedIndex.size(); ++i)
286 : {
287 : // get constrained rhs
288 : // modify block directly since set to zero afterwards
289 0 : block_type& val = rhs[constrainedIndex[i]];
290 0 : val *= frac;
291 :
292 : // split equally on all constraining indices
293 0 : for(size_t j=0; j < vConstrainingIndex.size(); ++j)
294 0 : rhs[vConstrainingIndex[j][i]] += val;
295 :
296 : // set rhs to zero for constrained index
297 0 : val = 0.0;
298 : }
299 0 : }
300 :
301 : template <typename TVector>
302 0 : void SplitAddRhs_OneSide(TVector& rhs,
303 : std::vector<size_t> & constrainedIndex,
304 : std::vector<std::vector<size_t> >& vConstrainingIndex)
305 : {
306 : typedef typename TVector::value_type block_type;
307 :
308 : // check number of indices passed
309 0 : for(size_t i = 0; i < vConstrainingIndex.size(); ++i)
310 : UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(),
311 : "Wrong number of indices.");
312 :
313 0 : for(size_t i = 0; i < constrainedIndex.size(); ++i)
314 : {
315 : // choose randomly the first dof to add whole
316 : // (must be the same as for row)
317 0 : const size_t addTo = vConstrainingIndex[0][i];
318 :
319 0 : block_type& val = rhs[constrainedIndex[i]];
320 0 : rhs[addTo] += val;
321 0 : val = 0.0;
322 : }
323 0 : }
324 :
325 :
326 : /**
327 : * @brief Extract algebra indices for constrained and constraining indices from DoF distribution
328 : *
329 : * @param dd DoF distribution
330 : * @param hgVrt the hanging vertex
331 : * @param vConstrainingVrt vector of constraining vertices
332 : * @param constrainedInd vector of algebra indices on hanging vertex
333 : * @param vConstrainingInd vector of (vector of constraining algebra indices) for constraining vertices
334 : */
335 0 : inline void get_algebra_indices(ConstSmartPtr<DoFDistribution> dd,
336 : ConstrainedVertex* hgVrt,
337 : std::vector<Vertex*>& vConstrainingVrt,
338 : std::vector<size_t>& constrainedInd,
339 : std::vector<std::vector<size_t> >& vConstrainingInd)
340 : {
341 : // get subset index
342 0 : const int si = dd->subset_handler()->get_subset_index(hgVrt);
343 :
344 : // get constraining vertices
345 0 : CollectConstraining(vConstrainingVrt, *dd->multi_grid(), hgVrt);
346 :
347 : // clear constrainedInd
348 : constrainedInd.clear();
349 :
350 : // resize constraining indices
351 : vConstrainingInd.clear();
352 0 : vConstrainingInd.resize(vConstrainingVrt.size());
353 :
354 0 : if (dd->grouped()) // block algebra (assuming constrainers have exactly the same functions as constrained)
355 : {
356 : // get indices for constrained vertex
357 0 : dd->inner_algebra_indices(hgVrt, constrainedInd, false);
358 :
359 : // get indices for constraining vertices
360 0 : for (size_t i = 0; i < vConstrainingVrt.size(); ++i)
361 0 : dd->inner_algebra_indices(vConstrainingVrt[i], vConstrainingInd[i], false);
362 : }
363 : else // scalar algebra
364 : {
365 : // get algebra indices for constrained and constraining vertices
366 0 : for (size_t fct = 0; fct < dd->num_fct(); ++fct)
367 : {
368 : // check that function is defined on subset
369 0 : if (!dd->is_def_in_subset(fct, si)) continue;
370 :
371 : // get indices for constrained vertex
372 0 : dd->inner_algebra_indices_for_fct(hgVrt, constrainedInd, false, fct);
373 :
374 : // get indices for constraining vertices
375 0 : for (size_t i = 0; i < vConstrainingVrt.size(); ++i)
376 : {
377 0 : const int siC = dd->subset_handler()->get_subset_index(vConstrainingVrt[i]);
378 :
379 : // check that function is defined on subset
380 0 : UG_COND_THROW(!dd->is_def_in_subset(fct, siC),
381 : "Function " << fct << " is defined for a constrained vertex, "
382 : "but not for one of its constraining vertices!");
383 :
384 0 : dd->inner_algebra_indices_for_fct(vConstrainingVrt[i], vConstrainingInd[i], false, fct);
385 : }
386 : }
387 : }
388 0 : }
389 :
390 : ////////////////////////////////////////////////////////////////////////////////
391 : ////////////////////////////////////////////////////////////////////////////////
392 : // Sym P1 Constraints
393 : ////////////////////////////////////////////////////////////////////////////////
394 : ////////////////////////////////////////////////////////////////////////////////
395 :
396 : template <typename TDomain, typename TAlgebra>
397 : void
398 0 : SymP1Constraints<TDomain,TAlgebra>::
399 : adjust_defect(vector_type& d, const vector_type& u,
400 : ConstSmartPtr<DoFDistribution> dd, int type, number time,
401 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
402 : const std::vector<number>* vScaleMass,
403 : const std::vector<number>* vScaleStiff)
404 : {
405 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
406 0 : UG_THROW("index-wise assemble routine is not "
407 : "implemented for SymP1Constraints \n");
408 :
409 : // storage for indices and vertices
410 : std::vector<std::vector<size_t> > vConstrainingInd;
411 : std::vector<size_t> constrainedInd;
412 : std::vector<Vertex*> vConstrainingVrt;
413 :
414 : // get begin end of hanging vertices
415 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
416 0 : iter = dd->begin<ConstrainedVertex>();
417 0 : iterEnd = dd->end<ConstrainedVertex>();
418 :
419 : // loop constrained vertices
420 0 : for(; iter != iterEnd; ++iter)
421 : {
422 : // get hanging vert
423 : ConstrainedVertex* hgVrt = *iter;
424 :
425 : // get algebra indices for constrained and constraining vertices
426 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
427 :
428 : // adapt rhs
429 0 : SplitAddRhs_Symmetric(d, constrainedInd, vConstrainingInd);
430 : }
431 0 : }
432 :
433 :
434 : template <typename TDomain, typename TAlgebra>
435 : void
436 0 : SymP1Constraints<TDomain,TAlgebra>::
437 : adjust_rhs(vector_type& rhs, const vector_type& u,
438 : ConstSmartPtr<DoFDistribution> dd, int type, number time)
439 : {
440 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
441 0 : UG_THROW("index-wise assemble routine is not "
442 : "implemented for SymP1Constraints \n");
443 :
444 : // storage for indices and vertices
445 : std::vector<std::vector<size_t> > vConstrainingInd;
446 : std::vector<size_t> constrainedInd;
447 : std::vector<Vertex*> vConstrainingVrt;
448 :
449 : // get begin end of hanging vertices
450 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
451 0 : iter = dd->begin<ConstrainedVertex>();
452 0 : iterEnd = dd->end<ConstrainedVertex>();
453 :
454 : // loop constrained vertices
455 0 : for(; iter != iterEnd; ++iter)
456 : {
457 : // get hanging vert
458 : ConstrainedVertex* hgVrt = *iter;
459 :
460 : // get algebra indices for constrained and constraining vertices
461 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
462 :
463 : // adapt rhs
464 0 : SplitAddRhs_Symmetric(rhs, constrainedInd, vConstrainingInd);
465 : }
466 0 : }
467 :
468 : template <typename TDomain, typename TAlgebra>
469 : void
470 0 : SymP1Constraints<TDomain,TAlgebra>::
471 : adjust_jacobian(matrix_type& J, const vector_type& u,
472 : ConstSmartPtr<DoFDistribution> dd, int type, number time,
473 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
474 : const number s_a0)
475 : {
476 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
477 0 : UG_THROW("index-wise assemble routine is not "
478 : "implemented for SymP1Constraints \n");
479 :
480 : // storage for indices and vertices
481 : std::vector<std::vector<size_t> > vConstrainingInd;
482 : std::vector<size_t> constrainedInd;
483 : std::vector<Vertex*> vConstrainingVrt;
484 :
485 : // get begin end of hanging vertices
486 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
487 0 : iter = dd->begin<ConstrainedVertex>();
488 0 : iterEnd = dd->end<ConstrainedVertex>();
489 :
490 : // loop constrained vertices
491 0 : for(; iter != iterEnd; ++iter)
492 : {
493 : // get hanging vert
494 : ConstrainedVertex* hgVrt = *iter;
495 :
496 : // get algebra indices for constrained and constraining vertices
497 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
498 :
499 : // Split using indices
500 0 : SplitAddRow_Symmetric(J, constrainedInd, vConstrainingInd);
501 :
502 : // set interpolation
503 0 : SetInterpolation(J, constrainedInd, vConstrainingInd, m_bAssembleLinearProblem);
504 : }
505 0 : }
506 :
507 : template <typename TDomain, typename TAlgebra>
508 : void
509 0 : SymP1Constraints<TDomain,TAlgebra>::
510 : adjust_linear(matrix_type& mat, vector_type& rhs,
511 : ConstSmartPtr<DoFDistribution> dd, int type, number time)
512 : {
513 0 : m_bAssembleLinearProblem = true;
514 :
515 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
516 0 : UG_THROW("index-wise assemble routine is not "
517 : "implemented for SymP1Constraints \n");
518 :
519 : // storage for indices and vertices
520 : std::vector<std::vector<size_t> > vConstrainingInd;
521 : std::vector<size_t> constrainedInd;
522 : std::vector<Vertex*> vConstrainingVrt;
523 :
524 : // get begin end of hanging vertices
525 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
526 0 : iter = dd->begin<ConstrainedVertex>();
527 0 : iterEnd = dd->end<ConstrainedVertex>();
528 :
529 : // loop constrained vertices
530 0 : for(; iter != iterEnd; ++iter)
531 : {
532 : // get hanging vert
533 : ConstrainedVertex* hgVrt = *iter;
534 :
535 : // get algebra indices for constrained and constraining vertices
536 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
537 :
538 : // Split using indices
539 0 : SplitAddRow_Symmetric(mat, constrainedInd, vConstrainingInd);
540 :
541 : // set interpolation
542 0 : SetInterpolation(mat, constrainedInd, vConstrainingInd, true);
543 :
544 : // adapt rhs
545 0 : SplitAddRhs_Symmetric(rhs, constrainedInd, vConstrainingInd);
546 : }
547 0 : }
548 :
549 : template <typename TDomain, typename TAlgebra>
550 : void
551 0 : SymP1Constraints<TDomain,TAlgebra>::
552 : adjust_solution(vector_type& u, ConstSmartPtr<DoFDistribution> dd,
553 : int type, number time)
554 : {
555 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
556 0 : UG_THROW("index-wise assemble routine is not "
557 : "implemented for SymP1Constraints \n");
558 :
559 : // storage for indices and vertices
560 : std::vector<std::vector<size_t> > vConstrainingInd;
561 : std::vector<size_t> constrainedInd;
562 : std::vector<Vertex*> vConstrainingVrt;
563 :
564 : // get begin end of hanging vertices
565 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
566 0 : iter = dd->begin<ConstrainedVertex>();
567 0 : iterEnd = dd->end<ConstrainedVertex>();
568 :
569 : // loop constraining edges
570 0 : for(; iter != iterEnd; ++iter)
571 : {
572 : // get hanging vert
573 : ConstrainedVertex* hgVrt = *iter;
574 :
575 : // get algebra indices for constrained and constraining vertices
576 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
577 :
578 : // Interpolate values
579 0 : InterpolateValues(u, constrainedInd, vConstrainingInd);
580 : }
581 0 : }
582 :
583 :
584 : template <typename TDomain, typename TAlgebra>
585 : void
586 0 : SymP1Constraints<TDomain,TAlgebra>::
587 : adjust_prolongation
588 : (
589 : matrix_type& P,
590 : ConstSmartPtr<DoFDistribution> ddFine,
591 : ConstSmartPtr<DoFDistribution> ddCoarse,
592 : int type,
593 : number time
594 : )
595 : {
596 0 : if (m_bAssembleLinearProblem) return;
597 :
598 0 : if (this->m_spAssTuner->single_index_assembling_enabled())
599 0 : UG_THROW("index-wise assemble routine is not "
600 : "implemented for SymP1Constraints \n");
601 :
602 : // storage for indices and vertices
603 : std::vector<std::vector<size_t> > vConstrainingInd;
604 : std::vector<size_t> constrainedInd;
605 : std::vector<Vertex*> vConstrainingVrt;
606 :
607 : // get begin end of hanging vertices
608 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
609 0 : iter = ddFine->begin<ConstrainedVertex>();
610 0 : iterEnd = ddFine->end<ConstrainedVertex>();
611 :
612 : // loop constrained vertices
613 0 : for(; iter != iterEnd; ++iter)
614 : {
615 : // get hanging vert
616 : ConstrainedVertex* hgVrt = *iter;
617 :
618 : // get algebra indices for constrained and constraining vertices
619 0 : get_algebra_indices(ddFine, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
620 :
621 : // set zero row
622 : size_t sz = constrainedInd.size();
623 0 : for (size_t i = 0; i < sz; ++i)
624 0 : SetRow(P, constrainedInd[i], 0.0);
625 : }
626 0 : }
627 :
628 :
629 : template <typename TDomain, typename TAlgebra>
630 : void
631 0 : SymP1Constraints<TDomain,TAlgebra>::
632 : adjust_restriction
633 : (
634 : matrix_type& R,
635 : ConstSmartPtr<DoFDistribution> ddCoarse,
636 : ConstSmartPtr<DoFDistribution> ddFine,
637 : int type,
638 : number time
639 : )
640 : {
641 :
642 0 : }
643 :
644 :
645 :
646 : template <typename TDomain, typename TAlgebra>
647 : void
648 0 : SymP1Constraints<TDomain,TAlgebra>::
649 : adjust_correction
650 : ( vector_type& u,
651 : ConstSmartPtr<DoFDistribution> dd,
652 : int type,
653 : number time
654 : )
655 : {
656 : //typedef typename vector_type::value_type block_type;
657 :
658 0 : if (this->m_spAssTuner->single_index_assembling_enabled())
659 0 : UG_THROW("index-wise assemble routine is not "
660 : "implemented for OneSideP1Constraints \n");
661 :
662 : // storage for indices and vertices
663 : std::vector<std::vector<size_t> > vConstrainingInd;
664 : std::vector<size_t> constrainedInd;
665 : std::vector<Vertex*> vConstrainingVrt;
666 :
667 : // get begin end of hanging vertices
668 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
669 0 : iter = dd->begin<ConstrainedVertex>();
670 0 : iterEnd = dd->end<ConstrainedVertex>();
671 :
672 : // loop constrained vertices
673 0 : for (; iter != iterEnd; ++iter)
674 : {
675 : // get hanging vert
676 : ConstrainedVertex* hgVrt = *iter;
677 :
678 : // get algebra indices for constrained and constraining vertices
679 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
680 :
681 : // set all entries corresponding to constrained dofs to zero
682 0 : for (size_t i = 0; i < constrainedInd.size(); ++i)
683 0 : u[constrainedInd[i]] = 0.0;
684 : }
685 0 : }
686 :
687 :
688 : ////////////////////////////////////////////////////////////////////////////////
689 : ////////////////////////////////////////////////////////////////////////////////
690 : // OneSide P1 Constraints
691 : ////////////////////////////////////////////////////////////////////////////////
692 : ////////////////////////////////////////////////////////////////////////////////
693 :
694 : template<int dim>
695 : struct SortVertexPos {
696 :
697 : SortVertexPos(SmartPtr<Domain<dim, MultiGrid, MGSubsetHandler> > spDomain)
698 : : m_aaPos(spDomain->position_accessor())
699 : {}
700 :
701 : inline bool operator() (Vertex* vrt1, Vertex* vrt2)
702 : {UG_THROW(dim <<" not implemented.");}
703 :
704 : protected:
705 : typename Domain<dim, MultiGrid, MGSubsetHandler>::position_accessor_type& m_aaPos;
706 : };
707 :
708 : template<>
709 : inline bool SortVertexPos<1>::operator() (Vertex* vrt1, Vertex* vrt2)
710 : {
711 : if(m_aaPos[vrt1][0] < m_aaPos[vrt2][0]) {
712 : return true;
713 : }
714 : return false;
715 : }
716 :
717 : template<>
718 : inline bool SortVertexPos<2>::operator() (Vertex* vrt1, Vertex* vrt2)
719 : {
720 : if(m_aaPos[vrt1][0] < m_aaPos[vrt2][0]) {
721 : return true;
722 : }
723 : else if(m_aaPos[vrt1][0] == m_aaPos[vrt2][0]) {
724 : if(m_aaPos[vrt1][1] < m_aaPos[vrt2][1])
725 : return true;
726 : }
727 : return false;
728 : }
729 :
730 : template<>
731 : inline bool SortVertexPos<3>::operator() (Vertex* vrt1, Vertex* vrt2)
732 : {
733 : if(m_aaPos[vrt1][0] < m_aaPos[vrt2][0]) {
734 : return true;
735 : }
736 : else if(m_aaPos[vrt1][0] == m_aaPos[vrt2][0]) {
737 : if(m_aaPos[vrt1][1] < m_aaPos[vrt2][1]){
738 : return true;
739 : }
740 : else if(m_aaPos[vrt1][1] == m_aaPos[vrt2][1]){
741 : if(m_aaPos[vrt1][2] < m_aaPos[vrt2][2])
742 : return true;
743 : }
744 : return false;
745 : }
746 : return false;
747 : }
748 :
749 :
750 : /**
751 : * @brief Extract DoF indices for constrained and constraining indices from DoF distribution
752 : *
753 : * One cannot simply use algebra indices as constrainers and constrained vertices
754 : * might not have the same number of functions defined on them. Mapping correct
755 : * indices is only possible through DoF indices.
756 : *
757 : * @param dd DoF distribution
758 : * @param hgVrt the hanging vertex
759 : * @param vConstrainingVrt vector of constraining vertices
760 : * @param constrainedInd vector of DoFs indices on hanging vertex
761 : * @param vConstrainingInd vector of (vector of constraining DoF indices) for constraining vertices
762 : * @param sortVertexPos sorting functional for constrainers
763 : */
764 : template <typename TDomain>
765 : inline void get_algebra_indices(ConstSmartPtr<DoFDistribution> dd,
766 : ConstrainedVertex* hgVrt,
767 : std::vector<Vertex*>& vConstrainingVrt,
768 : std::vector<size_t>& constrainedInd,
769 : std::vector<std::vector<size_t> >& vConstrainingInd,
770 : const SortVertexPos<TDomain::dim>& sortVertexPos)
771 : {
772 : // get subset index
773 : const int si = dd->subset_handler()->get_subset_index(hgVrt);
774 :
775 : // get constraining vertices
776 : CollectConstraining(vConstrainingVrt, *dd->multi_grid(), hgVrt);
777 :
778 : #ifdef UG_PARALLEL
779 : std::sort(vConstrainingVrt.begin(), vConstrainingVrt.end(), sortVertexPos);
780 : #endif
781 :
782 : // clear constrainedInd
783 : constrainedInd.clear();
784 :
785 : // resize constraining indices
786 : vConstrainingInd.clear();
787 : vConstrainingInd.resize(vConstrainingVrt.size());
788 :
789 : // get algebra indices for constrained and constraining vertices
790 : if (dd->grouped()) // block algebra (assuming constrainers have exactly the same functions as constrained)
791 : {
792 : // get indices for constrained vertex
793 : dd->inner_algebra_indices(hgVrt, constrainedInd, false);
794 :
795 : // get indices for constraining vertices
796 : for (size_t i = 0; i < vConstrainingVrt.size(); ++i)
797 : dd->inner_algebra_indices(vConstrainingVrt[i], vConstrainingInd[i], false);
798 : }
799 : else // scalar algebra
800 : {
801 : for (size_t fct = 0; fct < dd->num_fct(); fct++)
802 : {
803 : // check that function is defined on subset
804 : if (!dd->is_def_in_subset(fct, si)) continue;
805 :
806 : // get indices for constrained vertex
807 : dd->inner_algebra_indices_for_fct(hgVrt, constrainedInd, false, fct);
808 :
809 : // get indices for constraining vertices
810 : for (size_t i = 0; i < vConstrainingVrt.size(); ++i)
811 : {
812 : const int siC = dd->subset_handler()->get_subset_index(vConstrainingVrt[i]);
813 :
814 : // check that function is defined on subset
815 : UG_COND_THROW(!dd->is_def_in_subset(fct, siC),
816 : "Function " << fct << " is defined for a constrained vertex, "
817 : "but not for one of its constraining vertices!");
818 :
819 : dd->inner_algebra_indices_for_fct(vConstrainingVrt[i], vConstrainingInd[i], false, fct);
820 : }
821 : }
822 : }
823 : }
824 :
825 :
826 : template <typename TDomain, typename TAlgebra>
827 : void
828 0 : OneSideP1Constraints<TDomain,TAlgebra>::
829 : adjust_defect(vector_type& d, const vector_type& u,
830 : ConstSmartPtr<DoFDistribution> dd, int type, number time,
831 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
832 : const std::vector<number>* vScaleMass,
833 : const std::vector<number>* vScaleStiff)
834 : {
835 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
836 0 : UG_THROW("index-wise assemble routine is not "
837 : "implemented for OneSideP1Constraints \n");
838 :
839 : // storage for indices and vertices
840 : std::vector<std::vector<size_t> > vConstrainingInd;
841 : std::vector<size_t> constrainedInd;
842 : std::vector<Vertex*> vConstrainingVrt;
843 :
844 : #ifdef UG_PARALLEL
845 : SortVertexPos<TDomain::dim> sortVertexPos(this->approximation_space()->domain());
846 : #endif
847 :
848 : // get begin end of hanging vertices
849 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
850 0 : iter = dd->begin<ConstrainedVertex>();
851 0 : iterEnd = dd->end<ConstrainedVertex>();
852 :
853 : // loop constrained vertices
854 0 : for(; iter != iterEnd; ++iter)
855 : {
856 : // get hanging vert
857 : ConstrainedVertex* hgVrt = *iter;
858 :
859 : // get algebra indices for constrained and constraining vertices
860 : #ifdef UG_PARALLEL
861 : get_algebra_indices<TDomain>(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd, sortVertexPos);
862 : #else
863 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
864 : #endif
865 :
866 : // adapt rhs
867 0 : SplitAddRhs_OneSide(d, constrainedInd, vConstrainingInd);
868 : }
869 0 : }
870 :
871 :
872 : template <typename TDomain, typename TAlgebra>
873 : void
874 0 : OneSideP1Constraints<TDomain,TAlgebra>::
875 : adjust_rhs(vector_type& rhs, const vector_type& u,
876 : ConstSmartPtr<DoFDistribution> dd, int type, number time)
877 : {
878 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
879 0 : UG_THROW("index-wise assemble routine is not "
880 : "implemented for OneSideP1Constraints \n");
881 :
882 : // storage for indices and vertices
883 : std::vector<std::vector<size_t> > vConstrainingInd;
884 : std::vector<size_t> constrainedInd;
885 : std::vector<Vertex*> vConstrainingVrt;
886 :
887 : #ifdef UG_PARALLEL
888 : SortVertexPos<TDomain::dim> sortVertexPos(this->approximation_space()->domain());
889 : #endif
890 :
891 : // get begin end of hanging vertices
892 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
893 0 : iter = dd->begin<ConstrainedVertex>();
894 0 : iterEnd = dd->end<ConstrainedVertex>();
895 :
896 : // loop constrained vertices
897 0 : for(; iter != iterEnd; ++iter)
898 : {
899 : // get hanging vert
900 : ConstrainedVertex* hgVrt = *iter;
901 :
902 : // get algebra indices for constrained and constraining vertices
903 : #ifdef UG_PARALLEL
904 : get_algebra_indices<TDomain>(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd, sortVertexPos);
905 : #else
906 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
907 : #endif
908 :
909 : // adapt rhs
910 0 : SplitAddRhs_OneSide(rhs, constrainedInd, vConstrainingInd);
911 : }
912 0 : }
913 :
914 : template <typename TDomain, typename TAlgebra>
915 : void
916 0 : OneSideP1Constraints<TDomain,TAlgebra>::
917 : adjust_jacobian(matrix_type& J, const vector_type& u,
918 : ConstSmartPtr<DoFDistribution> dd, int type, number time,
919 : ConstSmartPtr<VectorTimeSeries<vector_type> > vSol,
920 : const number s_a0)
921 : {
922 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
923 0 : UG_THROW("index-wise assemble routine is not "
924 : "implemented for OneSideP1Constraints \n");
925 :
926 : // storage for indices and vertices
927 : std::vector<std::vector<size_t> > vConstrainingInd;
928 : std::vector<size_t> constrainedInd;
929 : std::vector<Vertex*> vConstrainingVrt;
930 :
931 : #ifdef UG_PARALLEL
932 : SortVertexPos<TDomain::dim> sortVertexPos(this->approximation_space()->domain());
933 : #endif
934 :
935 : // get begin end of hanging vertices
936 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
937 0 : iter = dd->begin<ConstrainedVertex>();
938 0 : iterEnd = dd->end<ConstrainedVertex>();
939 :
940 : // loop constrained vertices
941 0 : for(; iter != iterEnd; ++iter)
942 : {
943 : // get hanging vert
944 : ConstrainedVertex* hgVrt = *iter;
945 :
946 : // get algebra indices for constrained and constraining vertices
947 : #ifdef UG_PARALLEL
948 : get_algebra_indices<TDomain>(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd, sortVertexPos);
949 : #else
950 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
951 : #endif
952 :
953 : // Split using indices
954 0 : SplitAddRow_OneSide(J, constrainedInd, vConstrainingInd);
955 :
956 : // set interpolation
957 0 : SetInterpolation(J, constrainedInd, vConstrainingInd, m_bAssembleLinearProblem);
958 : }
959 0 : }
960 :
961 : template <typename TDomain, typename TAlgebra>
962 : void
963 0 : OneSideP1Constraints<TDomain,TAlgebra>::
964 : adjust_linear(matrix_type& mat, vector_type& rhs,
965 : ConstSmartPtr<DoFDistribution> dd, int type, number time)
966 : {
967 0 : m_bAssembleLinearProblem = true;
968 :
969 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
970 0 : UG_THROW("index-wise assemble routine is not "
971 : "implemented for OneSideP1Constraints \n");
972 :
973 : // storage for indices and vertices
974 : std::vector<std::vector<size_t> > vConstrainingInd;
975 : std::vector<size_t> constrainedInd;
976 : std::vector<Vertex*> vConstrainingVrt;
977 :
978 : #ifdef UG_PARALLEL
979 : SortVertexPos<TDomain::dim> sortVertexPos(this->approximation_space()->domain());
980 : #endif
981 :
982 : // get begin end of hanging vertices
983 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
984 0 : iter = dd->begin<ConstrainedVertex>();
985 0 : iterEnd = dd->end<ConstrainedVertex>();
986 :
987 : // loop constraining edges
988 0 : for(; iter != iterEnd; ++iter)
989 : {
990 : // get hanging vert
991 : ConstrainedVertex* hgVrt = *iter;
992 :
993 : // get algebra indices for constrained and constraining vertices
994 : #ifdef UG_PARALLEL
995 : get_algebra_indices<TDomain>(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd, sortVertexPos);
996 : #else
997 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
998 : #endif
999 :
1000 : // Split using indices
1001 0 : SplitAddRow_OneSide(mat, constrainedInd, vConstrainingInd);
1002 :
1003 : // Set interpolation
1004 0 : SetInterpolation(mat, constrainedInd, vConstrainingInd, true);
1005 :
1006 : // adapt rhs
1007 0 : SplitAddRhs_OneSide(rhs, constrainedInd, vConstrainingInd);
1008 : }
1009 0 : }
1010 :
1011 : template <typename TDomain, typename TAlgebra>
1012 : void
1013 0 : OneSideP1Constraints<TDomain,TAlgebra>::
1014 : adjust_solution(vector_type& u, ConstSmartPtr<DoFDistribution> dd,
1015 : int type, number time)
1016 : {
1017 0 : if(this->m_spAssTuner->single_index_assembling_enabled())
1018 0 : UG_THROW("index-wise assemble routine is not "
1019 : "implemented for OneSideP1Constraints \n");
1020 :
1021 : // storage for indices and vertices
1022 : std::vector<std::vector<size_t> > vConstrainingInd;
1023 : std::vector<size_t> constrainedInd;
1024 : std::vector<Vertex*> vConstrainingVrt;
1025 :
1026 : // get begin end of hanging vertices
1027 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
1028 0 : iter = dd->begin<ConstrainedVertex>();
1029 0 : iterEnd = dd->end<ConstrainedVertex>();
1030 :
1031 : // loop constraining edges
1032 0 : for(; iter != iterEnd; ++iter)
1033 : {
1034 : // get hanging vert
1035 : ConstrainedVertex* hgVrt = *iter;
1036 :
1037 : // get algebra indices for constrained and constraining vertices
1038 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
1039 :
1040 : // Interpolate values
1041 0 : InterpolateValues(u, constrainedInd, vConstrainingInd);
1042 : }
1043 0 : }
1044 :
1045 :
1046 :
1047 : template <typename TDomain, typename TAlgebra>
1048 : void
1049 0 : OneSideP1Constraints<TDomain,TAlgebra>::
1050 : adjust_prolongation
1051 : (
1052 : matrix_type& P,
1053 : ConstSmartPtr<DoFDistribution> ddFine,
1054 : ConstSmartPtr<DoFDistribution> ddCoarse,
1055 : int type,
1056 : number time
1057 : )
1058 : {
1059 0 : if (m_bAssembleLinearProblem) return;
1060 :
1061 0 : if (this->m_spAssTuner->single_index_assembling_enabled())
1062 0 : UG_THROW("index-wise assemble routine is not "
1063 : "implemented for SymP1Constraints \n");
1064 :
1065 : // storage for indices and vertices
1066 : std::vector<std::vector<size_t> > vConstrainingInd;
1067 : std::vector<size_t> constrainedInd;
1068 : std::vector<Vertex*> vConstrainingVrt;
1069 :
1070 : // get begin end of hanging vertices
1071 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
1072 0 : iter = ddFine->begin<ConstrainedVertex>();
1073 0 : iterEnd = ddFine->end<ConstrainedVertex>();
1074 :
1075 : // loop constrained vertices
1076 0 : for(; iter != iterEnd; ++iter)
1077 : {
1078 : // get hanging vert
1079 : ConstrainedVertex* hgVrt = *iter;
1080 :
1081 : // get algebra indices for constrained and constraining vertices
1082 0 : get_algebra_indices(ddFine, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
1083 :
1084 : // set zero row
1085 : size_t sz = constrainedInd.size();
1086 0 : for (size_t i = 0; i < sz; ++i)
1087 0 : SetRow(P, constrainedInd[i], 0.0);
1088 : }
1089 0 : }
1090 :
1091 :
1092 : template <typename TDomain, typename TAlgebra>
1093 : void
1094 0 : OneSideP1Constraints<TDomain,TAlgebra>::
1095 : adjust_restriction
1096 : (
1097 : matrix_type& R,
1098 : ConstSmartPtr<DoFDistribution> ddCoarse,
1099 : ConstSmartPtr<DoFDistribution> ddFine,
1100 : int type,
1101 : number time
1102 : )
1103 0 : {}
1104 :
1105 :
1106 :
1107 :
1108 : template <typename TDomain, typename TAlgebra>
1109 : void
1110 0 : OneSideP1Constraints<TDomain,TAlgebra>::
1111 : adjust_correction
1112 : ( vector_type& u,
1113 : ConstSmartPtr<DoFDistribution> dd,
1114 : int type,
1115 : number time
1116 : )
1117 : {
1118 : //typedef typename vector_type::value_type block_type;
1119 :
1120 0 : if (this->m_spAssTuner->single_index_assembling_enabled())
1121 0 : UG_THROW("index-wise assemble routine is not "
1122 : "implemented for OneSideP1Constraints \n");
1123 :
1124 : // storage for indices and vertices
1125 : std::vector<std::vector<size_t> > vConstrainingInd;
1126 : std::vector<size_t> constrainedInd;
1127 : std::vector<Vertex*> vConstrainingVrt;
1128 :
1129 : // get begin end of hanging vertices
1130 : DoFDistribution::traits<ConstrainedVertex>::const_iterator iter, iterEnd;
1131 0 : iter = dd->begin<ConstrainedVertex>();
1132 0 : iterEnd = dd->end<ConstrainedVertex>();
1133 :
1134 : // loop constrained vertices
1135 0 : for (; iter != iterEnd; ++iter)
1136 : {
1137 : // get hanging vert
1138 : ConstrainedVertex* hgVrt = *iter;
1139 :
1140 : // get algebra indices for constrained and constraining vertices
1141 0 : get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd);
1142 :
1143 : // set all entries corresponding to constrained dofs to zero
1144 0 : for (size_t i = 0; i < constrainedInd.size(); ++i)
1145 0 : u[constrainedInd[i]] = 0.0;
1146 : }
1147 0 : }
1148 :
1149 :
1150 : }; // namespace ug
1151 :
1152 :
1153 :
1154 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__CONSTRAINTS__CONTINUITY_CONSTRAINTS__P1_CONTINUITY_CONSTRAINTS_IMPL__ */
|
