Nektar++
LinearElasticSystem.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: LinearElasticSystem.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: LinearElasticSystem solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <algorithm>
36#include <iomanip>
37
41#include <LocalRegions/SegExp.h>
47
48#ifdef NEKTAR_USE_MPI
50#endif
51
52#ifdef NEKTAR_USE_PETSC
54#endif
55
56namespace Nektar
57{
58
61 "LinearElasticSystem", LinearElasticSystem::create);
62
63/*
64 * @brief Generate mapping between a simplex and the reference element.
65 *
66 * Mapping that requires the coordinates of the three vertices (x,y) of the
67 * triangular element and outputs the function coefficients that defines the
68 * mapping to the reference element for each coordinate (xfunc and yfunc)
69 *
70 * @param x co-ordinates of triangle/tetrahedron
71 * @param xf components of mapping
72 */
74{
75 int n = geom->GetNumVerts(), i, j;
76
77 DNekMat map(n, n, 1.0, eFULL);
78 DNekMat mapref(n, n, 1.0, eFULL);
79
80 // Extract coordinate information.
81 for (i = 0; i < n; ++i)
82 {
83 for (j = 0; j < n - 1; ++j)
84 {
85 map(j, i) = (*geom->GetVertex(i))[j];
86 }
87 }
88
89 // Set up reference triangle or tetrahedron mapping.
90 if (n == 3)
91 {
92 mapref(0, 0) = -1.0;
93 mapref(1, 0) = -1.0;
94 mapref(0, 1) = 1.0;
95 mapref(1, 1) = -1.0;
96 mapref(0, 2) = -1.0;
97 mapref(1, 2) = 1.0;
98 }
99 else if (n == 4)
100 {
101 mapref(0, 0) = -1.0;
102 mapref(1, 0) = -1.0;
103 mapref(2, 0) = -1.0;
104 mapref(0, 1) = 1.0;
105 mapref(1, 1) = -1.0;
106 mapref(2, 1) = -1.0;
107 mapref(0, 2) = -1.0;
108 mapref(1, 2) = 1.0;
109 mapref(2, 2) = -1.0;
110 mapref(0, 3) = -1.0;
111 mapref(1, 3) = -1.0;
112 mapref(2, 3) = 1.0;
113 }
114
115 map.Invert();
116
117 DNekMat newmap = mapref * map;
118 DNekMat mapred(n - 1, n - 1, 1.0, eFULL);
119
120 for (i = 0; i < n - 1; ++i)
121 {
122 for (j = 0; j < n - 1; ++j)
123 {
124 mapred(i, j) = newmap(i, j);
125 }
126 }
127
128 return mapred;
129}
130
131/**
132 * @brief Default constructor.
133 */
137 : EquationSystem(pSession, pGraph)
138{
139}
140
141/**
142 * @brief Set up the linear elasticity system.
143 *
144 * This routine loads the E and nu variables from the session file, creates a
145 * coupled assembly map and creates containers for the statically condensed
146 * block matrix system.
147 */
148void LinearElasticSystem::v_InitObject(bool DeclareFields)
149{
150 EquationSystem::v_InitObject(DeclareFields);
151
152 const int nVel = m_fields[0]->GetCoordim(0);
153 int n;
154
155 ASSERTL0(nVel > 1, "Linear elastic solver not set up for"
156 " this dimension (only 2D/3D supported).");
157
158 // Make sure that we have Young's modulus and Poisson ratio set.
159 m_session->LoadParameter("E", m_E, 1.00);
160 m_session->LoadParameter("nu", m_nu, 0.25);
161 m_session->LoadParameter("Beta", m_beta, 1.00);
162
163 // Create a coupled assembly map which allows us to tie u, v and w
164 // fields together.
166 std::dynamic_pointer_cast<MultiRegions::ContField>(m_fields[0]);
168 m_session, m_graph, u->GetLocalToGlobalMap(),
169 m_fields[0]->GetBndConditions(), m_fields);
170
171 // Figure out size of our new matrix systems by looping over all expansions
172 // and multiply number of coefficients by velocity components.
173 const int nEl = m_fields[0]->GetExpSize();
175
176 Array<OneD, unsigned int> sizeBnd(nEl);
177 Array<OneD, unsigned int> sizeInt(nEl);
178
179 // Allocate storage for boundary and interior map caches.
182
183 // Cache interior and boundary maps to avoid recalculating
184 // constantly. Should really be handled by manager in StdRegions but not
185 // really working yet.
186 for (n = 0; n < nEl; ++n)
187 {
188 exp = m_fields[0]->GetExp(n);
189 sizeBnd[n] = nVel * exp->NumBndryCoeffs();
190 sizeInt[n] = nVel * exp->GetNcoeffs() - sizeBnd[n];
191 exp->GetBoundaryMap(m_bmap[n]);
192 exp->GetInteriorMap(m_imap[n]);
193 }
194
195 // Create block matrix storage for the statically condensed system.
196 MatrixStorage blkmatStorage = eDIAGONAL;
198 sizeBnd, sizeBnd, blkmatStorage);
200 blkmatStorage);
202 blkmatStorage);
204 blkmatStorage);
205}
206
207/**
208 * @brief Build matrix system for linear elasticity equations.
209 *
210 * This routine constructs the matrix discretisation arising from the weak
211 * formulation of the linear elasticity equations. The resulting matrix system
212 * is then passed to LinearElasticSystem::SetStaticCondBlock in order to
213 * construct the statically condensed system.
214 *
215 * All of the matrices involved in the construction of the divergence of the
216 * stress tensor are Laplacian-like. We use the variable coefficient
217 * functionality within the library to multiply the Laplacian by the appropriate
218 * constants for all diagonal terms, and off-diagonal terms use
219 * LinearElasticSystem::BuildLaplacianIJMatrix to construct matrices containing
220 *
221 * \f[ \int \partial_{x_i} \phi_k \partial_{x_j} \phi_l \f]
222 *
223 * where mixed derivative terms are present. Symmetry (in terms of k,l) is
224 * exploited to avoid constructing this matrix repeatedly.
225 *
226 * @todo Make static condensation optional and construct full system instead.
227 */
229{
230 const int nEl = m_fields[0]->GetExpSize();
231 const int nVel = m_fields[0]->GetCoordim(0);
232
234 int n;
235
236 // Factors map for matrix keys.
238
239 // Calculate various constants
240 NekDouble mu = m_E * 0.5 / (1.0 + m_nu);
241 NekDouble lambda = m_E * m_nu / (1.0 + m_nu) / (1.0 - 2.0 * m_nu);
242
243 bool verbose = m_session->DefinesCmdLineArgument("verbose");
244 bool root = m_session->GetComm()->GetRank() == 0;
245
246 // Loop over each element and construct matrices.
247 if (nVel == 2)
248 {
249 for (n = 0; n < nEl; ++n)
250 {
251 exp = m_fields[0]->GetExp(n);
252 const int nPhys = exp->GetTotPoints();
253 Array<OneD, NekDouble> aArr(nPhys, lambda + 2 * mu);
254 Array<OneD, NekDouble> bArr(nPhys, mu);
255
256 StdRegions::VarCoeffMap varcoeffA, varcoeffD;
257 varcoeffA[StdRegions::eVarCoeffD00] = aArr;
258 varcoeffA[StdRegions::eVarCoeffD11] = bArr;
259 varcoeffD[StdRegions::eVarCoeffD00] = bArr;
260 varcoeffD[StdRegions::eVarCoeffD11] = aArr;
261
263 exp->DetShapeType(), *exp, factors,
264 varcoeffA);
266 exp->DetShapeType(), *exp, factors,
267 varcoeffD);
268
269 /*
270 * mat holds the linear operator [ A B ] acting on [ u ].
271 * [ C D ] [ v ]
272 */
274 mat[0][0] = exp->GenMatrix(matkeyA);
275 mat[0][1] = BuildLaplacianIJMatrix(1, 0, mu + lambda, exp);
276 mat[1][0] = mat[0][1];
277 mat[1][1] = exp->GenMatrix(matkeyD);
278
279 // Set up the statically condensed block for this element.
280 SetStaticCondBlock(n, exp, mat);
281
282 if (verbose && root)
283 {
284 std::cout << "\rBuilding matrix system: "
285 << (int)(100.0 * n / nEl) << "%" << std::flush;
286 }
287 }
288 }
289 else if (nVel == 3)
290 {
291 for (n = 0; n < nEl; ++n)
292 {
293 exp = m_fields[0]->GetExp(n);
294 const int nPhys = exp->GetTotPoints();
295 Array<OneD, NekDouble> aArr(nPhys, lambda + 2 * mu);
296 Array<OneD, NekDouble> bArr(nPhys, mu);
297
298 StdRegions::VarCoeffMap varcoeffA, varcoeffE, varcoeffI;
299 varcoeffA[StdRegions::eVarCoeffD00] = aArr;
300 varcoeffA[StdRegions::eVarCoeffD11] = bArr;
301 varcoeffA[StdRegions::eVarCoeffD22] = bArr;
302 varcoeffE[StdRegions::eVarCoeffD00] = bArr;
303 varcoeffE[StdRegions::eVarCoeffD11] = aArr;
304 varcoeffE[StdRegions::eVarCoeffD22] = bArr;
305 varcoeffI[StdRegions::eVarCoeffD00] = bArr;
306 varcoeffI[StdRegions::eVarCoeffD11] = bArr;
307 varcoeffI[StdRegions::eVarCoeffD22] = aArr;
308
310 exp->DetShapeType(), *exp, factors,
311 varcoeffA);
313 exp->DetShapeType(), *exp, factors,
314 varcoeffE);
316 exp->DetShapeType(), *exp, factors,
317 varcoeffI);
318
319 /*
320 * mat holds the linear operator [ A B C ] acting on [ u ].
321 * [ D E F ] [ v ]
322 * [ G H I ] [ w ]
323 */
325 mat[0][0] = exp->GenMatrix(matkeyA);
326 mat[0][1] = BuildLaplacianIJMatrix(1, 0, mu + lambda, exp);
327 mat[0][2] = BuildLaplacianIJMatrix(2, 0, mu + lambda, exp);
328
329 mat[1][0] = mat[0][1];
330 mat[1][1] = exp->GenMatrix(matkeyE);
331 mat[1][2] = BuildLaplacianIJMatrix(2, 1, mu + lambda, exp);
332
333 mat[2][0] = mat[0][2];
334 mat[2][1] = mat[1][2];
335 mat[2][2] = exp->GenMatrix(matkeyI);
336
337 // Set up the statically condensed block for this element.
338 SetStaticCondBlock(n, exp, mat);
339
340 if (verbose && root)
341 {
342 std::cout << "\rBuilding matrix system: "
343 << (int)(100.0 * n / nEl) << "%" << std::flush;
344 }
345 }
346 }
347
348 if (verbose && root)
349 {
350 std::cout << "\rBuilding matrix system: done." << std::endl;
351 }
352}
353
354/**
355 * @brief Generate summary at runtime.
356 */
358{
360
361 AddSummaryItem(s, "Young's modulus", m_E);
362 AddSummaryItem(s, "Poisson ratio", m_nu);
363}
364
365/**
366 * @brief Solve elliptic linear elastic system.
367 *
368 * The solve proceeds as follows:
369 *
370 * - Create a MultiRegions::GlobalLinSys object.
371 * - Evaluate a forcing function from the session file.
372 * - If the temperature term is enabled, evaluate this and add to forcing.
373 * - Apply Dirichlet boundary conditions.
374 * - Scatter forcing into the correct ordering according to the coupled
375 * assembly map.
376 * - Do the solve.
377 * - Scatter solution back to fields and backwards transform to physical space.
378 */
380{
381 int i, j, k, l, nv;
382 const int nVel = m_fields[0]->GetCoordim(0);
383
384 // Build initial matrix system.
386
387 // Now we've got the matrix system set up, create a GlobalLinSys
388 // object. We mask ourselves as LinearAdvectionReaction to create a full
389 // matrix instead of symmetric storage.
393
394 // Currently either direct or iterative static condensation is
395 // supported.
396 if (m_assemblyMap->GetGlobalSysSolnType() ==
398 {
402 }
403 else if (m_assemblyMap->GetGlobalSysSolnType() ==
405 {
410 }
411#ifdef NEKTAR_USE_PETSC
412 else if (m_assemblyMap->GetGlobalSysSolnType() ==
414 {
418 }
419#endif
420#ifdef NEKTAR_USE_MPI
421 else if (m_assemblyMap->GetGlobalSysSolnType() ==
423 {
427 }
428#endif
429
430 linSys->Initialise(m_assemblyMap);
431
432 const int nCoeffs = m_fields[0]->GetNcoeffs();
433
434 //
435 // -- Evaluate forcing functions
436 //
437
438 // Evaluate the forcing function from the XML file.
440 GetFunction("Forcing")->Evaluate(forcing);
441
442 // Add temperature term
443 std::string tempEval;
444 m_session->LoadSolverInfo("Temperature", tempEval, "None");
445
446 if (tempEval == "Jacobian")
447 {
448 // Allocate storage
450
451 for (nv = 0; nv < nVel; ++nv)
452 {
454 m_temperature[nv] =
456
457 for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
458 {
459 // Calculate element area
460 LocalRegions::ExpansionSharedPtr exp = m_fields[0]->GetExp(i);
461 LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
462 Array<OneD, NekDouble> jac = exp->GetMetricInfo()->GetJac(pkey);
463
464 int offset = m_fields[0]->GetPhys_Offset(i);
465
466 if (exp->GetMetricInfo()->GetGtype() ==
468 {
469 Vmath::Smul(exp->GetTotPoints(), m_beta, jac, 1,
470 tmp = m_temperature[nv] + offset, 1);
471 }
472 else
473 {
474 Vmath::Fill(exp->GetTotPoints(), m_beta * jac[0],
475 tmp = m_temperature[nv] + offset, 1);
476 }
477 }
478 m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
479 }
480 }
481 else if (tempEval == "Metric")
482 {
483 ASSERTL0((m_fields[0]->GetCoordim(0) == 2 &&
484 m_graph->GetAllQuadGeoms().size() == 0) ||
485 (m_fields[0]->GetCoordim(0) == 3 &&
486 m_graph->GetAllPrismGeoms().size() == 0 &&
487 m_graph->GetAllPyrGeoms().size() == 0 &&
488 m_graph->GetAllHexGeoms().size() == 0),
489 "LinearIdealMetric temperature only implemented for "
490 "two-dimensional triangular meshes or three-dimensional "
491 "tetrahedral meshes.");
492
495
496 for (nv = 0; nv < nVel; ++nv)
497 {
498 m_temperature[nv] =
501
502 for (i = 0; i < nVel; ++i)
503 {
504 m_stress[nv][i] =
506 }
507 }
508
509 for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
510 {
511 // Calculate element area
512 LocalRegions::ExpansionSharedPtr exp = m_fields[0]->GetExp(i);
513 LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
515 exp->GetMetricInfo()->GetDeriv(pkey);
516 int offset = m_fields[0]->GetPhys_Offset(i);
517
518 DNekMat i2rm = MappingIdealToRef(exp->GetGeom());
519
520 // Compute metric tensor
521 DNekMat jac(nVel, nVel, 0.0, eFULL);
522 DNekMat jacIdeal(nVel, nVel, 0.0, eFULL);
523 DNekMat metric(nVel, nVel, 0.0, eFULL);
524
525 for (j = 0; j < deriv[0][0].size(); ++j)
526 {
527 for (k = 0; k < nVel; ++k)
528 {
529 for (l = 0; l < nVel; ++l)
530 {
531 jac(l, k) = deriv[k][l][j];
532 }
533 }
534
535 jacIdeal = jac * i2rm;
536 metric = Transpose(jacIdeal) * jacIdeal;
537
538 // Compute eigenvalues/eigenvectors of metric tensor using
539 // ideal mapping.
540 char jobvl = 'N', jobvr = 'V';
541 int worklen = 8 * nVel, info;
542
543 DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
544 DNekMat evec(nVel, nVel, 0.0, eFULL);
545 DNekMat evecinv(nVel, nVel, 0.0, eFULL);
546 Array<OneD, NekDouble> vl(nVel * nVel);
547 Array<OneD, NekDouble> work(worklen);
548 Array<OneD, NekDouble> wi(nVel);
549
550 Lapack::Dgeev(jobvl, jobvr, nVel, metric.GetRawPtr(), nVel,
551 &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
552 &(evec.GetPtr())[0], nVel, &work[0], worklen,
553 info);
554
555 evecinv = evec;
556 evecinv.Invert();
557
558 // rescaling of the eigenvalues
559 for (nv = 0; nv < nVel; ++nv)
560 {
561 eval(nv, nv) = m_beta * (sqrt(eval(nv, nv)) - 1.0);
562 }
563
564 DNekMat beta = evec * eval * evecinv;
565
566 for (k = 0; k < nVel; ++k)
567 {
568 for (l = 0; l < nVel; ++l)
569 {
570 m_stress[k][l][offset + j] = beta(k, l);
571 }
572 }
573 }
574
575 if (deriv[0][0].size() != exp->GetTotPoints())
576 {
578 for (k = 0; k < nVel; ++k)
579 {
580 for (l = 0; l < nVel; ++l)
581 {
582 Vmath::Fill(exp->GetTotPoints(), m_stress[k][l][offset],
583 tmp = m_stress[k][l] + offset, 1);
584 }
585 }
586 }
587 }
588
589 // Calculate divergence of stress tensor.
591 for (i = 0; i < nVel; ++i)
592 {
593 for (j = 0; j < nVel; ++j)
594 {
595 m_fields[i]->PhysDeriv(j, m_stress[i][j], tmpderiv);
596 Vmath::Vadd(m_fields[i]->GetNpoints(), tmpderiv, 1,
597 m_temperature[i], 1, m_temperature[i], 1);
598 }
599
601 forcing[i], 1);
602 }
603 }
604 else if (tempEval != "None")
605 {
606 ASSERTL0(false, "Unknown temperature form: " + tempEval);
607 }
608
609 // Set up some temporary storage.
610 //
611 // - forCoeffs holds the forcing coefficients in a local ordering;
612 // however note that the ordering is different and dictated by the
613 // assembly map. We loop over each element, then the boundary degrees of
614 // freedom for u, boundary for v, followed by the interior for u and then
615 // interior for v.
616 // - inout holds the Dirichlet degrees of freedom in the local ordering,
617 // which have the boundary conditions imposed
618 Array<OneD, NekDouble> forCoeffs(nVel * nCoeffs, 0.0);
619 Array<OneD, NekDouble> inout(nVel * nCoeffs, 0.0);
620 Array<OneD, NekDouble> tmp(nCoeffs);
621
622 for (nv = 0; nv < nVel; ++nv)
623 {
624 // Take the inner product of the forcing function.
625 m_fields[nv]->IProductWRTBase(forcing[nv], tmp);
626
627 // Impose Dirichlet condition on field which should be initialised
628 Array<OneD, NekDouble> loc_inout = m_fields[nv]->UpdateCoeffs();
629 m_fields[nv]->ImposeDirichletConditions(loc_inout);
630
631 // Scatter forcing into RHS vector according to the ordering dictated in
632 // the comment above.
633 for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
634 {
635 int nBnd = m_bmap[i].size();
636 int nInt = m_imap[i].size();
637 int offset = m_fields[nv]->GetCoeff_Offset(i);
638
639 for (j = 0; j < nBnd; ++j)
640 {
641 forCoeffs[nVel * offset + nv * nBnd + j] =
642 -1.0 * tmp[offset + m_bmap[i][j]];
643 inout[nVel * offset + nv * nBnd + j] =
644 loc_inout[offset + m_bmap[i][j]];
645 }
646 for (j = 0; j < nInt; ++j)
647 {
648 forCoeffs[nVel * (offset + nBnd) + nv * nInt + j] =
649 -1.0 * tmp[offset + m_imap[i][j]];
650 }
651 }
652 }
653
654 //
655 // -- Perform solve
656 //
657 // Solve.
658 linSys->Solve(forCoeffs, inout, m_assemblyMap);
659
660 //
661 // -- Postprocess
662 //
663
664 // Scatter back to field degrees of freedom
665 for (nv = 0; nv < nVel; ++nv)
666 {
667 for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
668 {
669 int nBnd = m_bmap[i].size();
670 int nInt = m_imap[i].size();
671 int offset = m_fields[nv]->GetCoeff_Offset(i);
672
673 for (j = 0; j < nBnd; ++j)
674 {
675 m_fields[nv]->UpdateCoeffs()[offset + m_bmap[i][j]] =
676 inout[nVel * offset + nv * nBnd + j];
677 }
678 for (j = 0; j < nInt; ++j)
679 {
680 m_fields[nv]->UpdateCoeffs()[offset + m_imap[i][j]] =
681 inout[nVel * (offset + nBnd) + nv * nInt + j];
682 }
683 }
684 m_fields[nv]->BwdTrans(m_fields[nv]->GetCoeffs(),
685 m_fields[nv]->UpdatePhys());
686 }
687}
688
689/**
690 * @brief Given a block matrix for an element, construct its static condensation
691 * matrices.
692 *
693 * This routine essentially duplicates the logic present in many of the
694 * LocalRegions matrix generation routines to construct the statically condensed
695 * equivalent of mat to pass to the GlobalLinSys solver.
696 *
697 * @param n Element/block number.
698 * @param exp Pointer to expansion.
699 * @param mat Block matrix containing matrix operator.
700 */
702 const int n, const LocalRegions::ExpansionSharedPtr exp,
704{
705 int i, j, k, l;
706 const int nVel = mat.GetRows();
707 const int nB = exp->NumBndryCoeffs();
708 const int nI = exp->GetNcoeffs() - nB;
709 const int nBnd = exp->NumBndryCoeffs() * nVel;
710 const int nInt = exp->GetNcoeffs() * nVel - nBnd;
711 const MatrixStorage s = eFULL; // Maybe look into doing symmetric
712 // version of this?
713
722
723 for (i = 0; i < nVel; ++i)
724 {
725 for (j = 0; j < nVel; ++j)
726 {
727 // Boundary-boundary and boundary-interior
728 for (k = 0; k < nB; ++k)
729 {
730 for (l = 0; l < nB; ++l)
731 {
732 (*A)(k + i * nB, l + j * nB) =
733 (*mat[i][j])(m_bmap[n][k], m_bmap[n][l]);
734 }
735
736 for (l = 0; l < nI; ++l)
737 {
738 (*B)(k + i * nB, l + j * nI) =
739 (*mat[i][j])(m_bmap[n][k], m_imap[n][l]);
740 }
741 }
742
743 // Interior-boundary / interior-interior
744 for (k = 0; k < nI; ++k)
745 {
746 for (l = 0; l < nB; ++l)
747 {
748 (*C)(k + i * nI, l + j * nB) =
749 (*mat[i][j])(m_imap[n][k], m_bmap[n][l]);
750 }
751
752 for (l = 0; l < nI; ++l)
753 {
754 (*D)(k + i * nI, l + j * nI) =
755 (*mat[i][j])(m_imap[n][k], m_imap[n][l]);
756 }
757 }
758 }
759 }
760
761 // Construct static condensation matrices.
762 D->Invert();
763 (*B) = (*B) * (*D);
764 (*A) = (*A) - (*B) * (*C);
765
766 DNekScalMatSharedPtr tmp_mat;
767 m_schurCompl->SetBlock(
769 m_BinvD->SetBlock(
770 n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, B));
771 m_C->SetBlock(
772 n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, C));
773 m_Dinv->SetBlock(
774 n, n, tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, D));
775}
776
777/**
778 * @brief Construct a LaplacianIJ matrix for a given expansion.
779 *
780 * This routine constructs matrices whose entries contain the evaluation of
781 *
782 * \f[ \partial_{\tt k1} \phi_i \partial_{\tt k2} \phi_j \,dx \f]
783 */
785 const int k1, const int k2, const NekDouble scale,
787{
788 const int nCoeffs = exp->GetNcoeffs();
789 const int nPhys = exp->GetTotPoints();
790 int i;
791
792 DNekMatSharedPtr ret =
793 MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs, 0.0, eFULL);
794
795 Array<OneD, NekDouble> tmp2(nPhys);
796 Array<OneD, NekDouble> tmp3(nPhys);
797
798 for (i = 0; i < nCoeffs; ++i)
799 {
800 Array<OneD, NekDouble> tmp1(nCoeffs, 0.0);
801 tmp1[i] = 1.0;
802
803 exp->BwdTrans(tmp1, tmp2);
804 exp->PhysDeriv(k1, tmp2, tmp3);
805 exp->IProductWRTDerivBase(k2, tmp3, tmp1);
806
807 Vmath::Smul(nCoeffs, scale, &tmp1[0], 1,
808 &(ret->GetPtr())[0] + i * nCoeffs, 1);
809 }
810
811 return ret;
812}
813
815 std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
816 std::vector<std::string> &variables)
817{
818 const int nVel = m_fields[0]->GetCoordim(0);
819 const int nCoeffs = m_fields[0]->GetNcoeffs();
820
821 if (m_temperature.size() == 0)
822 {
823 return;
824 }
825
826 for (int i = 0; i < nVel; ++i)
827 {
828 Array<OneD, NekDouble> tFwd(nCoeffs);
829 m_fields[i]->FwdTrans(m_temperature[i], tFwd);
830 fieldcoeffs.push_back(tFwd);
831 variables.push_back("ThermStressDiv" + std::to_string(i));
832 }
833
834 if (m_stress.size() == 0)
835 {
836 return;
837 }
838
839 for (int i = 0; i < nVel; ++i)
840 {
841 for (int j = 0; j < nVel; ++j)
842 {
843 Array<OneD, NekDouble> tFwd(nCoeffs);
844 m_fields[i]->FwdTrans(m_stress[i][j], tFwd);
845 fieldcoeffs.push_back(tFwd);
846 variables.push_back("ThermStress" + std::to_string(i) +
847 std::to_string(j));
848 }
849 }
850}
851
852} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
NekDouble m_nu
Poisson ratio.
void v_InitObject(bool DeclareFields=true) override
Set up the linear elasticity system.
void BuildMatrixSystem()
Build matrix system for linear elasticity equations.
LinearElasticSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Default constructor.
void SetStaticCondBlock(const int n, const LocalRegions::ExpansionSharedPtr exp, Array< TwoD, DNekMatSharedPtr > &mat)
Given a block matrix for an element, construct its static condensation matrices.
Array< OneD, Array< OneD, NekDouble > > m_temperature
Storage for the temperature terms.
CoupledAssemblyMapSharedPtr m_assemblyMap
Assembly map for the coupled (u,v,w) system.
void v_DoSolve() override
Solve elliptic linear elastic system.
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables) override
NekDouble m_E
Young's modulus.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
Storage for the thermal stress terms.
static std::string className
Name of class.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
DNekMatSharedPtr BuildLaplacianIJMatrix(const int k1, const int k2, const NekDouble scale, LocalRegions::ExpansionSharedPtr exp)
Construct a LaplacianIJ matrix for a given expansion.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Generate summary at runtime.
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
A base class for describing how to solve specific equations.
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareFeld=true)
Initialisation object for EquationSystem.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
static void Dgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
Solve general real matrix eigenproblem.
Definition: Lapack.hpp:329
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
static PreconditionerSharedPtr NullPreconditionerSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:51
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
DNekMat MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285