MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
lscartesiancontrolpoint.cpp
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#include "lscartesiancontrolpoint.h"
8
9#include <fstream>
10#include <iomanip>
11#include "GEOM/geometry.h"
12#include "GEOM/geometryelement.h"
13#include "INCLUDE/maiaconstants.h"
14#include "MEMORY/alloc.h"
15#include "MEMORY/collector.h"
16#include "MEMORY/scratch.h"
17#include "UTIL/debug.h"
18#include "UTIL/timer.h"
19
20using namespace std;
21
22template <MInt nDim>
23LsControlPoint<nDim>::LsControlPoint() {
24 TRACE();
25 geometry = nullptr;
26}
27
28template <MInt nDim>
29LsControlPoint<nDim>::~LsControlPoint() {
30 TRACE();
31 geometry = nullptr;
32}
33
34template <MInt nDim>
35void LsControlPoint<nDim>::CtrlPnt1_Initialize(Geometry<nDim>* g, MInt NotUsed(dimensions)) {
36 TRACE();
37
38 geometry = g;
39 noElements = geometry->m_mbelements->size();
40 // initialize geometry at ref. geometry mid-point
41 for(MInt j = 0; j < nDim; j++)
42 CtrlPnt1_LagrangPnt[0][j] = 0.0;
43 // the other 3 points at bounding box maximum in the direction
44 for(MInt i = 0; i < 4; i++)
45 for(MInt j = 0; j < 3; j++)
46 CtrlPnt1_LagrangPnt[i][j] = 0.0;
47 for(MInt i = 1; i <= nDim; i++) {
48 for(MInt j = 0; j < nDim; j++) {
49 if(j == (i - 1)) CtrlPnt1_LagrangPnt[i][j] = geometry->m_mbminMax[j + nDim];
50 }
51 }
52 // Special pre-defined for 2D case
53 IF_CONSTEXPR(nDim == 2) CtrlPnt1_LagrangPnt[3][2] = 1.0;
54}
55
56template <MInt nDim>
57void LsControlPoint<nDim>::CtrlPnt1_InitPosition(MFloat* p) {
58 TRACE();
59
60 for(MInt i = 0; i <= 3; i++)
61 for(MInt j = 0; j < 3; j++) {
62 CtrlPnt1_LagrangPnt[i][j] += p[j];
63 }
64}
65
66template <MInt nDim>
67void LsControlPoint<nDim>::CtrlPnt1_InitOrientation(MFloat* u, MFloat* v, MFloat* w) {
68 TRACE();
69
70 // Transpose already
71 for(MInt i = 0; i < 3; i++) {
72 tr[0][i] = u[i];
73 tr[1][i] = v[i];
74 tr[2][i] = w[i];
75 }
76
77 MFloat** R = nullptr;
78 mAlloc(R, 3, 3, "R", F0, AT_);
79 if(!MATXINVT(tr, R)) {
80 for(MInt i = 0; i < 3; i++)
81 for(MInt j = 0; j < 3; j++)
82 R[i][j] = tr[i][j];
83 }
84 for(MInt i = 0; i < 3; i++)
85 for(MInt j = 0; j < 3; j++)
86 trinv[i][j] = R[i][j];
87 for(MInt i = 0; i < 4; i++) {
88 MFloat quvw[3] = {0, 0, 0}, qxyz[3] = {0, 0, 0};
89 for(MInt j = 0; j < 3; j++)
90 quvw[j] = CtrlPnt1_LagrangPnt[i][j];
91 for(MInt k = 0; k < 3; k++)
92 for(MInt j = 0; j < 3; j++) {
93 qxyz[k] += trinv[k][j] * quvw[j];
94 }
95 for(MInt j = 0; j < 3; j++)
96 CtrlPnt1_LagrangPnt[i][j] = qxyz[j];
97 }
98
99 mDeallocate(R);
100}
101
102template <MInt nDim>
103void LsControlPoint<nDim>::CtrlPnt1_UpdateTR() {
104 TRACE();
105
106 MFloat u[3], v[3], w[3];
107 for(MInt i = 0; i < 3; i++) {
108 u[i] = CtrlPnt1_LagrangPnt[1][i] - CtrlPnt1_LagrangPnt[0][i];
109 v[i] = CtrlPnt1_LagrangPnt[2][i] - CtrlPnt1_LagrangPnt[0][i];
110 w[i] = CtrlPnt1_LagrangPnt[3][i] - CtrlPnt1_LagrangPnt[0][i];
111 }
112 MFloat norm_u = norm(u);
113 MFloat norm_v = norm(v);
114 MFloat norm_w = norm(w);
115 for(MInt i = 0; i < 3; i++) {
116 u[i] = u[i] / norm_u;
117 v[i] = v[i] / norm_v;
118 w[i] = w[i] / norm_w;
119 }
120 for(MInt j = 0; j < 3; j++) {
121 tr[0][j] = u[j];
122 tr[1][j] = v[j];
123 tr[2][j] = w[j];
124 }
125}
126
127// Warning vector q will be changed
128template <MInt nDim>
129void LsControlPoint<nDim>::CtrlPnt1_quvw(MFloat* q) {
130 TRACE();
131
132 CtrlPnt1_UpdateTR();
133 MFloat quvw[3] = {0, 0, 0};
134 for(MInt i = 0; i < 3; i++)
135 q[i] = q[i] - CtrlPnt1_LagrangPnt[0][i];
136 for(MInt i = 0; i < 3; i++)
137 for(MInt j = 0; j < 3; j++) {
138 quvw[i] += tr[i][j] * q[j];
139 }
140 for(MInt i = 0; i < 3; i++)
141 q[i] = quvw[i];
142}
143
144template <MInt nDim>
145void LsControlPoint<nDim>::CtrlPnt1_CtrlPntToSTL(const MChar* FILENAME) {
146 TRACE();
147
148 IF_CONSTEXPR(nDim == 2) {
149 cerr << "CtrlPntToSTL is not valid for 2D." << endl;
150 return;
151 }
152 // Update TR
153 CtrlPnt1_UpdateTR();
154 // Find invert TR
155 MFloat** R = nullptr;
156 mAlloc(R, 3, 3, "R", F0, AT_);
157
158 if(!MATXINVT(tr, R)) {
159 for(MInt i = 0; i < 3; i++)
160 for(MInt j = 0; j < 3; j++)
161 R[i][j] = tr[i][j];
162 }
163 for(MInt i = 0; i < 3; i++)
164 for(MInt j = 0; j < 3; j++)
165 trinv[i][j] = R[i][j];
166
167 ofstream fo(FILENAME);
168 fo << "solid CTRL-PNT" << endl;
169 fo << scientific;
170 for(MInt e = 0; e < noElements; e++) {
171 fo << setw(15) << "facet normal";
172 MFloat nuvw[3] = {0, 0, 0}, nxyz[3] = {0, 0, 0};
173 for(MInt i = 0; i < 3; i++) {
174 nuvw[i] = geometry->mbelements[e].m_normal[i];
175 }
176 for(MInt i = 0; i < 3; i++)
177 for(MInt j = 0; j < 3; j++) {
178 nxyz[i] += trinv[i][j] * nuvw[j];
179 }
180 for(MInt i = 0; i < 3; i++) {
181 fo << setw(20) << nxyz[i] << '\t';
182 }
183 fo << endl;
184 fo << setw(15) << "outer loop" << endl;
185 for(MInt v = 0; v < 3; v++) {
186 MFloat quvw[3] = {0, 0, 0}, qxyz[3] = {0, 0, 0};
187 fo << setw(15) << "vertex";
188 for(MInt i = 0; i < 3; i++)
189 quvw[i] = geometry->mbelements[e].m_vertices[v][i];
190 for(MInt i = 0; i < 3; i++)
191 for(MInt j = 0; j < 3; j++) {
192 qxyz[i] += trinv[i][j] * quvw[j];
193 }
194 for(MInt i = 0; i < 3; i++)
195 qxyz[i] += CtrlPnt1_LagrangPnt[0][i];
196 for(MInt i = 0; i < 3; i++) {
197 fo << setw(20) << qxyz[i];
198 }
199 fo << endl;
200 }
201 fo << setw(15) << "endloop" << endl;
202 fo << setw(15) << "endfacet" << endl;
203 }
204 fo << "endsolid" << endl;
205 fo.close();
206
207 mDeallocate(R);
208}
209
210template <MInt nDim>
211void LsControlPoint<nDim>::CtrlPnt2_Initialize(Geometry<nDim>* g, MInt NotUsed(dimensions)) {
212 TRACE();
213
214 geometry = g;
215 noElements = geometry->m_mbelements->size();
216}
217
218template <MInt nDim>
219void LsControlPoint<nDim>::CtrlPnt2_InitPosition(MFloat* p) {
220 TRACE();
221
222 geometry->MoveAllMBElementVertex(p);
223}
224
225template <MInt nDim>
226void LsControlPoint<nDim>::CtrlPnt2_InitOrientation(MFloat* u, MFloat* v, MFloat* w) {
227 TRACE();
228
229 // Transpose already
230 for(MInt i = 0; i < 3; i++) {
231 tr[0][i] = u[i];
232 tr[1][i] = v[i];
233 tr[2][i] = w[i];
234 }
235 MFloat** R = nullptr;
236 mAlloc(R, 3, 3, "R", F0, AT_);
237
238 if(!MATXINVT(tr, R)) {
239 for(MInt i = 0; i < 3; i++)
240 for(MInt j = 0; j < 3; j++)
241 R[i][j] = tr[i][j];
242 }
243 for(MInt i = 0; i < 3; i++)
244 for(MInt j = 0; j < 3; j++)
245 trinv[i][j] = R[i][j];
246 mDeallocate(R);
247 for(MInt e = 0; e < noElements; e++) {
248 for(MInt vt = 0; vt < nDim; vt++) {
249 MFloat quvw[3] = {0, 0, 0}, qxyz[3] = {0, 0, 0};
250 for(MInt i = 0; i < nDim; i++)
251 quvw[i] = geometry->mbelements[e].m_vertices[vt][i];
252 for(MInt i = 0; i < 3; i++)
253 for(MInt j = 0; j < 3; j++) {
254 qxyz[i] += trinv[i][j] * quvw[j];
255 }
256 geometry->ReplaceMBElementVertex(e, vt, qxyz);
257 }
258 }
259 CtrlPnt2_UpdateAllNormalVector();
260}
261
262template <MInt nDim>
263void LsControlPoint<nDim>::CtrlPnt2_getMBElementVertex(MInt e, MInt v, MFloat* r) {
264 TRACE();
265
266 for(MInt i = 0; i < nDim; i++)
267 r[i] = geometry->mbelements[e].m_vertices[v][i];
268}
269
270template <MInt nDim>
271void LsControlPoint<nDim>::CtrlPnt2_MoveMBElementVertex(MInt e, MInt v, MFloat* dx) {
272 TRACE();
273
274 geometry->MoveMBElementVertex(e, v, dx);
275}
276
277template <MInt nDim>
278void LsControlPoint<nDim>::CtrlPnt2_Update() {
279 TRACE();
280
281 geometry->UpdateMBBoundingBox();
282 geometry->UpdateADT();
283}
284
285template <MInt nDim>
286void LsControlPoint<nDim>::CtrlPnt2_UpdateAllNormalVector() {
287 TRACE();
288
289 for(MInt e = 0; e < noElements; e++)
290 geometry->UpdateMBNormalVector(e);
291}
292
293template <MInt nDim>
294void LsControlPoint<nDim>::CtrlPnt2_UpdateNormalVector(MInt e) {
295 TRACE();
296
297 geometry->UpdateMBNormalVector(e);
298}
299
300template <MInt nDim>
301void LsControlPoint<nDim>::CtrlPnt2_InitPosition(MFloat* p, MInt bcId) {
302 TRACE();
303
304 for(MInt e = 0; e < noElements; e++)
305 if(geometry->mbelements[e].m_bndCndId == bcId)
306 for(MInt v = 0; v < nDim; v++)
307 geometry->MoveMBElementVertex(e, v, p);
308}
309
310template <MInt nDim>
311void LsControlPoint<nDim>::CtrlPnt2_InitOrientation(MFloat* u, MFloat* v, MFloat* w, MInt bcId) {
312 TRACE();
313
314 // Transpose already
315 for(MInt i = 0; i < 3; i++) {
316 tr[0][i] = u[i];
317 tr[1][i] = v[i];
318 tr[2][i] = w[i];
319 }
320 MFloat** R = nullptr;
321 mAlloc(R, 3, 3, "R", F0, AT_);
322
323 if(!MATXINVT(tr, R)) {
324 for(MInt i = 0; i < 3; i++)
325 for(MInt j = 0; j < 3; j++)
326 R[i][j] = tr[i][j];
327 }
328 for(MInt i = 0; i < 3; i++)
329 for(MInt j = 0; j < 3; j++)
330 trinv[i][j] = R[i][j];
331 // for(MInt i=0;i<3;i++) {delete [] R[i];} delete [] R; R=nullptr;
332 mDeallocate(R);
333 for(MInt e = 0; e < noElements; e++) {
334 if(geometry->mbelements[e].m_bndCndId == bcId) {
335 for(MInt vt = 0; vt < nDim; vt++) {
336 MFloat quvw[3] = {0, 0, 0}, qxyz[3] = {0, 0, 0};
337 for(MInt i = 0; i < nDim; i++)
338 quvw[i] = geometry->mbelements[e].m_vertices[vt][i];
339 for(MInt i = 0; i < 3; i++)
340 for(MInt j = 0; j < 3; j++) {
341 qxyz[i] += trinv[i][j] * quvw[j];
342 }
343 geometry->ReplaceMBElementVertex(e, vt, qxyz);
344 }
345 }
346 }
347 CtrlPnt2_UpdateAllNormalVector();
348}
349
350template <MInt nDim>
351void LsControlPoint<nDim>::CtrlPnt2_MoveMBElementVertex(MInt e, MInt v, MFloat* dx, MInt bcId) {
352 TRACE();
353
354 if(geometry->mbelements[e].m_bndCndId == bcId) geometry->MoveMBElementVertex(e, v, dx);
355}
356
357template <MInt nDim>
358void LsControlPoint<nDim>::CtrlPnt2_shiftSTL(MInt bcId, MFloat* dx) {
359 TRACE();
360
361 // Movable STL
362 for(MInt e = 0; e < CtrlPnt2_getNoElements(); e++) {
363 for(MInt v = 0; v < nDim; v++) {
364 // Move control points to the new time
365 CtrlPnt2_MoveMBElementVertex(e, v, dx, bcId);
366 }
367 }
368 CtrlPnt2_Update();
369}
370
371template <MInt nDim>
372void LsControlPoint<nDim>::CtrlPnt2_rotateSTL(MInt bcId, MFloat* dphi, MFloat* ori) {
373 TRACE();
374
375 IF_CONSTEXPR(nDim != 3) {
376 mTerm(1, AT_, "Untested in 2D!");
377 // nPos[0] = ori[0] - cos(dphi[0]) * (oPos[0] - ori[0]) + sin(dphi[0]) * (oPos[0] - ori[0]);
378 // nPos[1] = ori[1] - sin(dphi[0]) * (oPos[0] - ori[0]) - cos(dphi[0]) * (oPos[0] - ori[0]);
379 }
380
381 MFloat oPos[nDim];
382 MFloat nPos[nDim];
383 MFloat dx[nDim];
384
385 MFloat rotMat[3][3];
386 rotMat[0][0] = (cos(dphi[1]) * cos(dphi[0]));
387 rotMat[0][1] = (cos(dphi[1]) * sin(dphi[0]));
388 rotMat[0][2] = (-sin(dphi[1]));
389 rotMat[1][0] = (sin(dphi[2]) * sin(dphi[1]) * cos(dphi[0]) - cos(dphi[2]) * sin(dphi[0]));
390 rotMat[1][1] = (sin(dphi[2]) * sin(dphi[1]) * sin(dphi[0]) + cos(dphi[2]) * cos(dphi[0]));
391 rotMat[1][2] = (sin(dphi[2]) * cos(dphi[1]));
392 rotMat[2][0] = (cos(dphi[2]) * sin(dphi[1]) * cos(dphi[0]) + sin(dphi[2]) * sin(dphi[0]));
393 rotMat[2][1] = (cos(dphi[2]) * sin(dphi[1]) * sin(dphi[0]) - sin(dphi[2]) * cos(dphi[0]));
394 rotMat[2][2] = (cos(dphi[2]) * cos(dphi[1]));
395
396 for(MInt e = 0; e < CtrlPnt2_getNoElements(); e++) {
397 if(geometry->mbelements[e].m_bndCndId != bcId) continue;
398
399 for(MInt v = 0; v < nDim; v++) {
400 for(MInt n = 0; n < nDim; n++) {
401 oPos[n] = geometry->mbelements[e].m_vertices[v][n];
402 }
403
404 nPos[0] = ori[0] + rotMat[0][0] * (oPos[0] - ori[0]) + rotMat[0][1] * (oPos[1] - ori[1])
405 + rotMat[0][2] * (oPos[2] - ori[2]);
406 nPos[1] = ori[1] + rotMat[1][0] * (oPos[0] - ori[0]) + rotMat[1][1] * (oPos[1] - ori[1])
407 + rotMat[1][2] * (oPos[2] - ori[2]);
408 nPos[2] = ori[2] + rotMat[2][0] * (oPos[0] - ori[0]) + rotMat[2][1] * (oPos[1] - ori[1])
409 + rotMat[2][2] * (oPos[2] - ori[2]);
410
411 for(MInt n = 0; n < nDim; n++) {
412 dx[n] = nPos[n] - oPos[n];
413 }
414
415 // Move control points to the new time
416 CtrlPnt2_MoveMBElementVertex(e, v, dx, bcId);
417 }
418 geometry->UpdateMBNormalVector(e);
419 }
420 CtrlPnt2_Update();
421}
422
423template <MInt nDim>
424MBool LsControlPoint<nDim>::MATXINVT(MFloat T[3][3], MFloat** R) {
425 TRACE();
426
427 MFloat** a = nullptr;
428 MFloat** y = nullptr;
429 MFloat col[3];
430 MInt indx[3];
431 MFloat d;
432
433 mAlloc(a, 3, 3, "a", F0, AT_);
434 mAlloc(y, 3, 3, "y", F0, AT_);
435
436 for(MInt i = 0; i < 3; i++)
437 for(MInt j = 0; j < 3; j++)
438 a[i][j] = T[i][j];
439 if(!LUDECMP(a, indx, d)) {
440 return false;
441 }
442 for(MInt j = 0; j < 3; j++) {
443 for(MInt i = 0; i < 3; i++)
444 col[i] = 0.0;
445 col[j] = 1.0;
446 LUBKSB(a, indx, col);
447 for(MInt i = 0; i < 3; i++)
448 y[i][j] = col[i];
449 }
450 for(MInt i = 0; i < 3; i++)
451 for(MInt j = 0; j < 3; j++)
452 R[i][j] = y[i][j];
453 mDeallocate(a);
454 mDeallocate(y);
455
456 return true;
457}
458
459template <MInt nDim>
460MBool LsControlPoint<nDim>::LUDECMP(MFloat** a, MInt* indx, MFloat& d) {
461 TRACE();
462
463 const MFloat TINY = 1.0e-20;
464 MInt imax = std::numeric_limits<MInt>::min();
465 // matrix size
466 static const MInt n = 3;
467 // scale factor
468 MFloat vv[3];
469 for(MInt i = 0; i < n; i++) {
470 MFloat big = 0.0;
471 for(MInt j = 0; j < n; j++) {
472 MFloat temp = fabs(a[i][j]);
473 if(temp > big) big = temp;
474 }
475 if(approx(big, 0.0, MFloatEps)) {
476 return false;
477 }
478 vv[i] = 1.0 / big;
479 }
480 for(MInt j = 0; j < n; j++) {
481 for(MInt i = 0; i < j; i++) {
482 MFloat sum = a[i][j];
483 for(MInt k = 0; k < i; k++)
484 sum -= a[i][k] * a[k][j];
485 a[i][j] = sum;
486 }
487 MFloat big = 0.0;
488 for(MInt i = j; i < n; i++) {
489 MFloat sum = a[i][j];
490 for(MInt k = 0; k < j; k++)
491 sum -= a[i][k] * a[k][j];
492 a[i][j] = sum;
493 MFloat dum = vv[i] * fabs(sum);
494 if(dum >= big) {
495 big = dum;
496 imax = i;
497 }
498 }
499 if(j != imax) {
500 for(MInt k = 0; k < n; k++) {
501 MFloat dum = a[imax][k];
502 a[imax][k] = a[j][k];
503 a[j][k] = dum;
504 }
505 d = -d;
506 vv[imax] = vv[j];
507 }
508 indx[j] = imax;
509 if(approx(a[j][j], 0.0, MFloatEps)) a[j][j] = TINY;
510 if(j != n - 1) {
511 MFloat dum = 1.0 / (a[j][j]);
512 for(MInt i = j + 1; i < n; i++)
513 a[i][j] *= dum;
514 }
515 }
516
517 return true;
518}
519
520template <MInt nDim>
521void LsControlPoint<nDim>::LUBKSB(MFloat** a, MInt* indx, MFloat* b) {
522 TRACE();
523
524 MInt i, ii = 0, ip, j;
525 MFloat sum;
526 MInt n = 3;
527 for(i = 0; i < n; i++) {
528 ip = indx[i];
529 sum = b[ip];
530 b[ip] = b[i];
531 if(ii != 0)
532 for(j = ii - 1; j < i; j++)
533 sum -= a[i][j] * b[j];
534 else if(!approx(sum, 0.0, MFloatEps))
535 ii = i + 1;
536 b[i] = sum;
537 }
538 for(i = n - 1; i >= 0; i--) {
539 sum = b[i];
540 for(j = i + 1; j < n; j++)
541 sum -= a[i][j] * b[j];
542 b[i] = sum / a[i][i];
543 }
544}
545
546template <MInt nDim>
547void LsControlPoint<nDim>::CtrlPnt2_CtrlPntToSTL(const MChar* FILENAME, MInt bcId) {
548 TRACE();
549
550 IF_CONSTEXPR(nDim == 2) {
551 MFloat middle[3] = {0.0, 0.0, 0.0};
552 ofstream fo(FILENAME);
553 MInt noElem = 0;
554 // compute middle
555 for(MInt i = 0; i < 2; i++) {
556 middle[i] = F0;
557 }
558 for(MInt e = 0; e < noElements; e++) {
559 if(geometry->mbelements[e].m_bndCndId == bcId) {
560 noElem++;
561 for(MInt v = 0; v < 2; v++) {
562 for(MInt i = 0; i < 2; i++) {
563 middle[i] += geometry->mbelements[e].m_vertices[0][i];
564 }
565 }
566 }
567 }
568 for(MInt i = 0; i < 2; i++) {
569 middle[i] = middle[i] / (noElem * 2);
570 }
571
572 fo << "solid CTRL-PNT" << endl;
573 fo << scientific;
574 for(MInt e = 0; e < noElements; e++) {
575 if(geometry->mbelements[e].m_bndCndId == bcId) {
576 fo << setw(15) << "facet normal";
577 for(MInt i = 0; i < 2; i++) {
578 fo << setw(20) << 0.0;
579 }
580 fo << setw(20) << 1.0;
581 fo << endl;
582 fo << setw(15) << "outer loop" << endl;
583 for(MInt v = 0; v < 2; v++) {
584 fo << setw(15) << "vertex";
585 for(MInt i = 0; i < 2; i++) {
586 fo << setw(20) << geometry->mbelements[e].m_vertices[v][i];
587 }
588 fo << setw(20) << F0;
589 fo << endl;
590 }
591 fo << setw(15) << "vertex";
592 for(MInt i = 0; i < 3; i++) {
593 fo << setw(20) << middle[i];
594 }
595 fo << endl;
596 fo << setw(15) << "endloop" << endl;
597 fo << setw(15) << "endfacet" << endl;
598 }
599 }
600 fo << "endsolid" << endl;
601 fo.close();
602 }
603 else {
604 ofstream fo(FILENAME);
605 fo << "solid CTRL-PNT" << endl;
606 fo << scientific;
607 for(MInt e = 0; e < noElements; e++) {
608 if(geometry->mbelements[e].m_bndCndId == bcId) {
609 fo << setw(15) << "facet normal";
610 for(MInt i = 0; i < 3; i++) {
611 fo << setw(20) << geometry->mbelements[e].m_normal[i];
612 }
613 fo << endl;
614 fo << setw(15) << "outer loop" << endl;
615 for(MInt v = 0; v < 3; v++) {
616 fo << setw(15) << "vertex";
617 for(MInt i = 0; i < 3; i++) {
618 fo << setw(20) << geometry->mbelements[e].m_vertices[v][i];
619 }
620 fo << endl;
621 }
622 fo << setw(15) << "endloop" << endl;
623 fo << setw(15) << "endfacet" << endl;
624 }
625 }
626 fo << "endsolid" << endl;
627 fo.close();
628 }
629}
630
631template <MInt nDim>
632void LsControlPoint<nDim>::CtrlPnt2_CtrlPntToGNUPlot(const MChar* FILENAME) {
633 TRACE();
634
635 ofstream fo(FILENAME);
636 for(MInt e = 0; e < noElements; e++) {
637 for(MInt v = 0; v < nDim; v++) {
638 for(MInt i = 0; i < nDim; i++) {
639 fo << geometry->mbelements[e].m_vertices[v][i] << '\t';
640 }
641 fo << endl;
642 }
643 }
644 fo.close();
645}
646
647template <MInt nDim>
648MFloat LsControlPoint<nDim>::norm(MFloat* r) {
649 return sqrt(pow(r[0], 2) + pow(r[1], 2) + pow(r[2], 2));
650}
651
652// Explicit instantiations for 2D and 3D
653template class LsControlPoint<2>;
654template class LsControlPoint<3>;
mAlloc
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
Definition: alloc.h:173
mDeallocate
MBool mDeallocate(T *&a)
deallocates the memory previously allocated for element 'a'
Definition: alloc.h:544
Geometry::m_mbelements
Collector< element< nDim > > * m_mbelements
Definition: geometry.h:217
LsControlPoint::CtrlPnt2_Initialize
void CtrlPnt2_Initialize(Geometry< nDim > *g, MInt dimensions)
Definition: lscartesiancontrolpoint.cpp:211
LsControlPoint::CtrlPnt2_UpdateNormalVector
void CtrlPnt2_UpdateNormalVector(MInt e)
Definition: lscartesiancontrolpoint.cpp:294
LsControlPoint::CtrlPnt1_CtrlPntToSTL
void CtrlPnt1_CtrlPntToSTL(const MChar *FILNAME)
Definition: lscartesiancontrolpoint.cpp:145
LsControlPoint::LUBKSB
void LUBKSB(MFloat **a, MInt *indx, MFloat *b)
Definition: lscartesiancontrolpoint.cpp:521
LsControlPoint::CtrlPnt2_InitOrientation
void CtrlPnt2_InitOrientation(MFloat *u, MFloat *v, MFloat *w)
Definition: lscartesiancontrolpoint.cpp:226
LsControlPoint::CtrlPnt2_CtrlPntToSTL
void CtrlPnt2_CtrlPntToSTL(const MChar *FILENAME, MInt bcId)
Definition: lscartesiancontrolpoint.cpp:547
LsControlPoint::CtrlPnt2_getMBElementVertex
void CtrlPnt2_getMBElementVertex(MInt e, MInt v, MFloat *r)
Definition: lscartesiancontrolpoint.cpp:263
LsControlPoint::CtrlPnt2_CtrlPntToGNUPlot
void CtrlPnt2_CtrlPntToGNUPlot(const MChar *FILENAME)
Definition: lscartesiancontrolpoint.cpp:632
LsControlPoint::CtrlPnt1_InitOrientation
void CtrlPnt1_InitOrientation(MFloat *u, MFloat *v, MFloat *w)
Definition: lscartesiancontrolpoint.cpp:67
LsControlPoint::CtrlPnt1_quvw
void CtrlPnt1_quvw(MFloat *q)
Definition: lscartesiancontrolpoint.cpp:129
LsControlPoint::CtrlPnt1_InitPosition
void CtrlPnt1_InitPosition(MFloat *p)
Definition: lscartesiancontrolpoint.cpp:57
LsControlPoint::CtrlPnt2_shiftSTL
void CtrlPnt2_shiftSTL(MInt bcId, MFloat *dx)
Definition: lscartesiancontrolpoint.cpp:358
LsControlPoint::norm
MFloat norm(MFloat *r)
Definition: lscartesiancontrolpoint.cpp:648
LsControlPoint::CtrlPnt2_rotateSTL
void CtrlPnt2_rotateSTL(MInt bcId, MFloat *dphi, MFloat *ori)
Definition: lscartesiancontrolpoint.cpp:372
LsControlPoint::CtrlPnt2_InitPosition
void CtrlPnt2_InitPosition(MFloat *p)
Definition: lscartesiancontrolpoint.cpp:219
LsControlPoint::CtrlPnt2_MoveMBElementVertex
void CtrlPnt2_MoveMBElementVertex(MInt e, MInt v, MFloat *dx)
Definition: lscartesiancontrolpoint.cpp:271
LsControlPoint::LUDECMP
MBool LUDECMP(MFloat **a, MInt *indx, MFloat &d)
Definition: lscartesiancontrolpoint.cpp:460
LsControlPoint::MATXINVT
MBool MATXINVT(MFloat T[3][3], MFloat **R)
Definition: lscartesiancontrolpoint.cpp:424
LsControlPoint::CtrlPnt1_Initialize
void CtrlPnt1_Initialize(Geometry< nDim > *g, MInt dimensions)
Definition: lscartesiancontrolpoint.cpp:35
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
approx
MBool approx(const T &, const U &, const T)
Definition: functions.h:272
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
MChar
char MChar
Definition: maiatypes.h:56
a
Definition: contexttypes.h:19
y
define array structures
Definition: geometrycontexttypes.h:14