Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: LinearElasticSystem solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <algorithm>
37 #include <iomanip>
38 
39 #include <LocalRegions/MatrixKey.h>
40 #include <LocalRegions/SegExp.h>
50 
51 #ifdef NEKTAR_USE_MPI
53 #endif
54 
55 #ifdef NEKTAR_USE_PETSC
57 #endif
58 
59 using namespace std;
60 
61 namespace Nektar
62 {
63 
64 string LinearElasticSystem::className = GetEquationSystemFactory().
65  RegisterCreatorFunction("LinearElasticSystem",
66  LinearElasticSystem::create);
67 
68 /*
69  * @brief Generate mapping between a simplex and the reference element.
70  *
71  * Mapping that requires the coordinates of the three vertices (x,y) of the
72  * triangular element and outputs the function coefficients that defines the
73  * mapping to the reference element for each coordinate (xfunc and yfunc)
74  *
75  * @param x co-ordinates of triangle/tetrahedron
76  * @param xf components of mapping
77  */
79 {
80  int n = geom->GetNumVerts(), i, j;
81 
82  DNekMat map (n, n, 1.0, eFULL);
83  DNekMat mapref(n, n, 1.0, eFULL);
84 
85  // Extract coordinate information.
86  for (i = 0; i < n; ++i)
87  {
88  for (j = 0; j < n-1; ++j)
89  {
90  map(j,i) = (*geom->GetVertex(i))[j];
91  }
92  }
93 
94  // Set up reference triangle or tetrahedron mapping.
95  if (n == 3)
96  {
97  mapref(0,0) = -1.0; mapref(1,0) = -1.0;
98  mapref(0,1) = 1.0; mapref(1,1) = -1.0;
99  mapref(0,2) = -1.0; mapref(1,2) = 1.0;
100  }
101  else if (n == 4)
102  {
103  mapref(0,0) = -1.0; mapref(1,0) = -1.0; mapref(2,0) = -1.0;
104  mapref(0,1) = 1.0; mapref(1,1) = -1.0; mapref(2,1) = -1.0;
105  mapref(0,2) = -1.0; mapref(1,2) = 1.0; mapref(2,2) = -1.0;
106  mapref(0,3) = -1.0; mapref(1,3) = -1.0; mapref(2,3) = 1.0;
107  }
108 
109  map.Invert();
110 
111  DNekMat newmap = mapref * map;
112  DNekMat mapred(n-1, n-1, 1.0, eFULL);
113 
114  for (i = 0; i < n-1; ++i)
115  {
116  for (j = 0; j < n-1; ++j)
117  {
118  mapred(i,j) = newmap(i,j);
119  }
120  }
121 
122  return mapred;
123 }
124 
125 
126 /**
127  * @brief Default constructor.
128  */
129 LinearElasticSystem::LinearElasticSystem(
130  const LibUtilities::SessionReaderSharedPtr& pSession)
131  : EquationSystem(pSession)
132 {
133 }
134 
135 /**
136  * @brief Set up the linear elasticity system.
137  *
138  * This routine loads the E and nu variables from the session file, creates a
139  * coupled assembly map and creates containers for the statically condensed
140  * block matrix system.
141  */
143 {
145 
146  const int nVel = m_fields[0]->GetCoordim(0);
147  int n;
148 
149  ASSERTL0(nVel > 1, "Linear elastic solver not set up for"
150  " this dimension (only 2D/3D supported).");
151 
152  // Make sure that we have Young's modulus and Poisson ratio set.
153  m_session->LoadParameter("E", m_E, 1.00);
154  m_session->LoadParameter("nu", m_nu, 0.25);
155  m_session->LoadParameter("Beta", m_beta, 1.00);
156 
157  // Create a coupled assembly map which allows us to tie u, v and w
158  // fields together.
159  if (nVel == 2)
160  {
161  MultiRegions::ContField2DSharedPtr u = boost::dynamic_pointer_cast<
165  m_graph,
166  u->GetLocalToGlobalMap(),
167  m_fields[0]->GetBndConditions(),
168  m_fields);
169  }
170 
171  if (nVel == 3)
172  {
173  MultiRegions::ContField3DSharedPtr u = boost::dynamic_pointer_cast<
177  m_graph,
178  u->GetLocalToGlobalMap(),
179  m_fields[0]->GetBndConditions(),
180  m_fields);
181  }
182 
183  // Figure out size of our new matrix systems by looping over all expansions
184  // and multiply number of coefficients by velocity components.
185  const int nEl = m_fields[0]->GetExpSize();
187 
188  Array<OneD, unsigned int> sizeBnd(nEl);
189  Array<OneD, unsigned int> sizeInt(nEl);
190 
191  // Allocate storage for boundary and interior map caches.
194 
195  // Cache interior and boundary maps to avoid recalculating
196  // constantly. Should really be handled by manager in StdRegions but not
197  // really working yet.
198  for (n = 0; n < nEl; ++n)
199  {
200  exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
201  sizeBnd[n] = nVel * exp->NumBndryCoeffs();
202  sizeInt[n] = nVel * exp->GetNcoeffs() - sizeBnd[n];
203  exp->GetBoundaryMap(m_bmap[n]);
204  exp->GetInteriorMap(m_imap[n]);
205  }
206 
207  // Create block matrix storage for the statically condensed system.
208  MatrixStorage blkmatStorage = eDIAGONAL;
210  sizeBnd, sizeBnd, blkmatStorage);
212  sizeBnd, sizeInt, blkmatStorage);
214  sizeInt, sizeBnd, blkmatStorage);
216  sizeInt, sizeInt, blkmatStorage);
217 }
218 
219 /**
220  * @brief Build matrix system for linear elasticity equations.
221  *
222  * This routine constructs the matrix discretisation arising from the weak
223  * formulation of the linear elasticity equations. The resulting matrix system
224  * is then passed to LinearElasticSystem::SetStaticCondBlock in order to
225  * construct the statically condensed system.
226  *
227  * All of the matrices involved in the construction of the divergence of the
228  * stress tensor are Laplacian-like. We use the variable coefficient
229  * functionality within the library to multiply the Laplacian by the appropriate
230  * constants for all diagonal terms, and off-diagonal terms use
231  * LinearElasticSystem::BuildLaplacianIJMatrix to construct matrices containing
232  *
233  * \f[ \int \partial_{x_i} \phi_k \partial_{x_j} \phi_l \f]
234  *
235  * where mixed derivative terms are present. Symmetry (in terms of k,l) is
236  * exploited to avoid constructing this matrix repeatedly.
237  *
238  * @todo Make static condensation optional and construct full system instead.
239  */
241 {
242  const int nEl = m_fields[0]->GetExpSize();
243  const int nVel = m_fields[0]->GetCoordim(0);
244 
246  int n;
247 
248  // Factors map for matrix keys.
250 
251  // Calculate various constants
252  NekDouble mu = m_E * 0.5 / (1.0 + m_nu);
253  NekDouble lambda = m_E * m_nu / (1.0 + m_nu) / (1.0 - 2.0*m_nu);
254 
255  bool verbose = m_session->DefinesCmdLineArgument("verbose");
256  bool root = m_session->GetComm()->GetRank() == 0;
257 
258  // Loop over each element and construct matrices.
259  if (nVel == 2)
260  {
261  for (n = 0; n < nEl; ++n)
262  {
263  exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
264  const int nPhys = exp->GetTotPoints();
265  Array<OneD, NekDouble> aArr(nPhys, lambda + 2*mu);
266  Array<OneD, NekDouble> bArr(nPhys, mu);
267 
268  StdRegions::VarCoeffMap varcoeffA, varcoeffD;
269  varcoeffA[StdRegions::eVarCoeffD00] = aArr;
270  varcoeffA[StdRegions::eVarCoeffD11] = bArr;
271  varcoeffD[StdRegions::eVarCoeffD00] = bArr;
272  varcoeffD[StdRegions::eVarCoeffD11] = aArr;
273 
275  exp->DetShapeType(),
276  *exp, factors, varcoeffA);
278  exp->DetShapeType(),
279  *exp, factors, varcoeffD);
280 
281  /*
282  * mat holds the linear operator [ A B ] acting on [ u ].
283  * [ C D ] [ v ]
284  */
286  mat[0][0] = exp->GenMatrix(matkeyA);
287  mat[0][1] = BuildLaplacianIJMatrix(1, 0, mu+lambda, exp);
288  mat[1][0] = mat[0][1];
289  mat[1][1] = exp->GenMatrix(matkeyD);
290 
291  // Set up the statically condensed block for this element.
292  SetStaticCondBlock(n, exp, mat);
293 
294  if (verbose && root)
295  {
296  cout << "\rBuilding matrix system: "
297  << (int)(100.0 * n / nEl) << "%" << flush;
298  }
299  }
300  }
301  else if (nVel == 3)
302  {
303  for (n = 0; n < nEl; ++n)
304  {
305  exp = m_fields[0]->GetExp(m_fields[0]->GetOffset_Elmt_Id(n));
306  const int nPhys = exp->GetTotPoints();
307  Array<OneD, NekDouble> aArr(nPhys, lambda + 2*mu);
308  Array<OneD, NekDouble> bArr(nPhys, mu);
309 
310  StdRegions::VarCoeffMap varcoeffA, varcoeffE, varcoeffI;
311  varcoeffA[StdRegions::eVarCoeffD00] = aArr;
312  varcoeffA[StdRegions::eVarCoeffD11] = bArr;
313  varcoeffA[StdRegions::eVarCoeffD22] = bArr;
314  varcoeffE[StdRegions::eVarCoeffD00] = bArr;
315  varcoeffE[StdRegions::eVarCoeffD11] = aArr;
316  varcoeffE[StdRegions::eVarCoeffD22] = bArr;
317  varcoeffI[StdRegions::eVarCoeffD00] = bArr;
318  varcoeffI[StdRegions::eVarCoeffD11] = bArr;
319  varcoeffI[StdRegions::eVarCoeffD22] = aArr;
320 
322  exp->DetShapeType(),
323  *exp, factors, varcoeffA);
325  exp->DetShapeType(),
326  *exp, factors, varcoeffE);
328  exp->DetShapeType(),
329  *exp, factors, varcoeffI);
330 
331  /*
332  * mat holds the linear operator [ A B C ] acting on [ u ].
333  * [ D E F ] [ v ]
334  * [ G H I ] [ w ]
335  */
337  mat[0][0] = exp->GenMatrix(matkeyA);
338  mat[0][1] = BuildLaplacianIJMatrix(1, 0, mu + lambda, exp);
339  mat[0][2] = BuildLaplacianIJMatrix(2, 0, mu + lambda, exp);
340 
341  mat[1][0] = mat[0][1];
342  mat[1][1] = exp->GenMatrix(matkeyE);
343  mat[1][2] = BuildLaplacianIJMatrix(2, 1, mu + lambda, exp);
344 
345  mat[2][0] = mat[0][2];
346  mat[2][1] = mat[1][2];
347  mat[2][2] = exp->GenMatrix(matkeyI);
348 
349  // Set up the statically condensed block for this element.
350  SetStaticCondBlock(n, exp, mat);
351 
352  if (verbose && root)
353  {
354  cout << "\rBuilding matrix system: "
355  << (int)(100.0 * n / nEl) << "%" << flush;
356  }
357  }
358  }
359 
360  if (verbose && root)
361  {
362  cout << "\rBuilding matrix system: done." << endl;
363  }
364 }
365 
366 /**
367  * @brief Generate summary at runtime.
368  */
370 {
372 
373  AddSummaryItem(s, "Young's modulus", m_E);
374  AddSummaryItem(s, "Poisson ratio", m_nu);
375 }
376 
377 /**
378  * @brief Solve elliptic linear elastic system.
379  *
380  * The solve proceeds as follows:
381  *
382  * - Create a MultiRegions::GlobalLinSys object.
383  * - Evaluate a forcing function from the session file.
384  * - If the temperature term is enabled, evaluate this and add to forcing.
385  * - Apply Dirichlet boundary conditions.
386  * - Scatter forcing into the correct ordering according to the coupled
387  * assembly map.
388  * - Do the solve.
389  * - Scatter solution back to fields and backwards transform to physical space.
390  */
392 {
393 
394  int i, j, k, l, nv;
395  const int nVel = m_fields[0]->GetCoordim(0);
396 
397  // Build initial matrix system.
399 
400  // Now we've got the matrix system set up, create a GlobalLinSys
401  // object. We mask ourselves as LinearAdvectionReaction to create a full
402  // matrix instead of symmetric storage.
406 
407  // Currently either direct or iterative static condensation is
408  // supported.
409  if (m_assemblyMap->GetGlobalSysSolnType() == MultiRegions::eDirectStaticCond)
410  {
411  linSys = MemoryManager<
412  MultiRegions::GlobalLinSysDirectStaticCond>::AllocateSharedPtr(
413  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
414  m_assemblyMap);
415  }
416  else if (m_assemblyMap->GetGlobalSysSolnType() ==
418  {
419  linSys = MemoryManager<
421  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
423  }
424 #ifdef NEKTAR_USE_PETSC
425  else if (m_assemblyMap->GetGlobalSysSolnType() ==
427  {
428  linSys = MemoryManager<
429  MultiRegions::GlobalLinSysPETScStaticCond>::AllocateSharedPtr(
430  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
431  m_assemblyMap);
432  }
433 #endif
434 #ifdef NEKTAR_USE_MPI
435  else if (m_assemblyMap->GetGlobalSysSolnType() ==
437  {
438  linSys = MemoryManager<
439  MultiRegions::GlobalLinSysXxtStaticCond>::AllocateSharedPtr(
440  key, m_fields[0], m_schurCompl, m_BinvD, m_C, m_Dinv,
441  m_assemblyMap);
442  }
443 #endif
444 
445  linSys->Initialise(m_assemblyMap);
446 
447  const int nCoeffs = m_fields[0]->GetNcoeffs();
448  const int nGlobDofs = m_assemblyMap->GetNumGlobalCoeffs();
449 
450  //
451  // -- Evaluate forcing functions
452  //
453 
454  // Evaluate the forcing function from the XML file.
455  Array<OneD, Array<OneD, NekDouble> > forcing(nVel);
456  EvaluateFunction(forcing, "Forcing");
457 
458  // Add temperature term
459  string tempEval;
460  m_session->LoadSolverInfo("Temperature", tempEval, "None");
461 
462  if (tempEval == "Jacobian")
463  {
464  // Allocate storage
466 
467  for (nv = 0; nv < nVel; ++nv)
468  {
471  m_fields[nv]->GetNpoints());
472 
473  for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
474  {
475  // Calculate element area
477  m_fields[0]->GetExp(i);
479  exp->GetPointsKeys();
481  exp->GetMetricInfo()->GetJac(pkey);
482 
483  int offset = m_fields[0]->GetPhys_Offset(i);
484 
485  if (exp->GetMetricInfo()->GetGtype() ==
487  {
488  Vmath::Smul(exp->GetTotPoints(), m_beta, jac, 1,
489  tmp = m_temperature[nv] + offset, 1);
490  }
491  else
492  {
493  Vmath::Fill(exp->GetTotPoints(), m_beta*jac[0],
494  tmp = m_temperature[nv] + offset, 1);
495  }
496  }
497  m_fields[nv]->PhysDeriv(nv, m_temperature[nv], forcing[nv]);
498  }
499  }
500  else if (tempEval == "Metric")
501  {
502  ASSERTL0((m_fields[0]->GetCoordim(0) == 2 &&
503  m_graph->GetAllQuadGeoms().size() == 0) ||
504  (m_fields[0]->GetCoordim(0) == 3 &&
505  m_graph->GetAllPrismGeoms().size() == 0 &&
506  m_graph->GetAllPyrGeoms ().size() == 0 &&
507  m_graph->GetAllHexGeoms ().size() == 0),
508  "LinearIdealMetric temperature only implemented for "
509  "two-dimensional triangular meshes or three-dimensional "
510  "tetrahedral meshes.");
511 
514 
515  for (nv = 0; nv < nVel; ++nv)
516  {
518  m_fields[nv]->GetNpoints(), 0.0);
519  m_stress[nv] = Array<OneD, Array<OneD, NekDouble> >(nVel);
520 
521  for (i = 0; i < nVel; ++i)
522  {
523  m_stress[nv][i] = Array<OneD, NekDouble>(
524  m_fields[nv]->GetNpoints(), 0.0);
525  }
526  }
527 
528  for (i = 0; i < m_fields[0]->GetExpSize(); ++i)
529  {
530  // Calculate element area
532  m_fields[0]->GetExp(i);
533  LibUtilities::PointsKeyVector pkey = exp->GetPointsKeys();
535  exp->GetMetricInfo()->GetDeriv(pkey);
536  int offset = m_fields[0]->GetPhys_Offset(i);
537 
538  DNekMat i2rm = MappingIdealToRef(exp->GetGeom());
539 
540  // Compute metric tensor
541  DNekMat jac (nVel, nVel, 0.0, eFULL);
542  DNekMat jacIdeal(nVel, nVel, 0.0, eFULL);
543  DNekMat metric (nVel, nVel, 0.0, eFULL);
544 
545  for (j = 0; j < deriv[0][0].num_elements(); ++j)
546  {
547  for (k = 0; k < nVel; ++k)
548  {
549  for (l = 0; l < nVel; ++l)
550  {
551  jac(l,k) = deriv[k][l][j];
552  }
553  }
554 
555  jacIdeal = jac * i2rm;
556  metric = Transpose(jacIdeal) * jacIdeal;
557 
558  // Compute eigenvalues/eigenvectors of metric tensor using
559  // ideal mapping.
560  char jobvl = 'N', jobvr = 'V';
561  int worklen = 8*nVel, info;
562 
563  DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
564  DNekMat evec (nVel, nVel, 0.0, eFULL);
565  DNekMat evecinv(nVel, nVel, 0.0, eFULL);
566  Array<OneD, NekDouble> vl (nVel*nVel);
567  Array<OneD, NekDouble> work(worklen);
568  Array<OneD, NekDouble> wi (nVel);
569 
570  Lapack::Dgeev(jobvl, jobvr, nVel, metric.GetRawPtr(), nVel,
571  &(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
572  &(evec.GetPtr())[0], nVel,
573  &work[0], worklen, info);
574 
575  evecinv = evec;
576  evecinv.Invert();
577 
578  // rescaling of the eigenvalues
579  for (nv = 0; nv < nVel; ++nv)
580  {
581  eval(nv,nv) = m_beta * (sqrt(eval(nv,nv)) - 1.0);
582  }
583 
584  DNekMat beta = evec * eval * evecinv;
585 
586  for (k = 0; k < nVel; ++k)
587  {
588  for (l = 0; l < nVel; ++l)
589  {
590  m_stress[k][l][offset+j] = beta(k,l);
591  }
592  }
593  }
594 
595  if (deriv[0][0].num_elements() != exp->GetTotPoints())
596  {
598  for (k = 0; k < nVel; ++k)
599  {
600  for (l = 0; l < nVel; ++l)
601  {
602  Vmath::Fill(
603  exp->GetTotPoints(), m_stress[k][l][offset],
604  tmp = m_stress[k][l] + offset, 1);
605  }
606  }
607  }
608  }
609 
610  // Calculate divergence of stress tensor.
612  for (i = 0; i < nVel; ++i)
613  {
614  for (j = 0; j < nVel; ++j)
615  {
616  m_fields[i]->PhysDeriv(j, m_stress[i][j], tmpderiv);
617  Vmath::Vadd (m_fields[i]->GetNpoints(), tmpderiv, 1,
618  m_temperature[i], 1, m_temperature[i], 1);
619  }
620 
622  m_temperature[i], 1, forcing[i], 1);
623  }
624  }
625  else if (tempEval != "None")
626  {
627  ASSERTL0(false, "Unknown temperature form: " + tempEval);
628  }
629 
630  // Set up some temporary storage.
631  //
632  // - forCoeffs holds the forcing coefficients in a local ordering;
633  // however note that the ordering is different and dictated by the
634  // assembly map. We loop over each element, then the boundary degrees of
635  // freedom for u, boundary for v, followed by the interior for u and then
636  // interior for v.
637  // - rhs is the global assembly of forCoeffs.
638  // - inout holds the Dirichlet degrees of freedom in the global ordering,
639  // which have been assembled from the boundary expansion.
640  Array<OneD, NekDouble> forCoeffs(nVel * nCoeffs, 0.0);
641  Array<OneD, NekDouble> inout (nGlobDofs, 0.0);
642  Array<OneD, NekDouble> rhs (nGlobDofs, 0.0);
643 
644  for (nv = 0; nv < nVel; ++nv)
645  {
646  // Take the inner product of the forcing function.
647  Array<OneD, NekDouble> tmp(nCoeffs);
648  m_fields[nv]->IProductWRTBase_IterPerExp(forcing[nv], tmp);
649 
650  // Scatter forcing into RHS vector according to the ordering dictated in
651  // the comment above.
652  for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
653  {
654  int nBnd = m_bmap[i].num_elements();
655  int nInt = m_imap[i].num_elements();
656  int offset = m_fields[nv]->GetCoeff_Offset(i);
657 
658  for (j = 0; j < nBnd; ++j)
659  {
660  forCoeffs[nVel*offset + nv*nBnd + j] =
661  tmp[offset+m_bmap[i][j]];
662  }
663  for (j = 0; j < nInt; ++j)
664  {
665  forCoeffs[nVel*(offset + nBnd) + nv*nInt + j] =
666  tmp[offset+m_imap[i][j]];
667  }
668  }
669  }
670 
671  // -- Impose Dirichlet boundary conditions.
672 
673  // First try to do parallel assembly: the intention here is that Dirichlet
674  // values at some edges/vertices need to be communicated to processes which
675  // don't contain the entire boundary region. See
676  // ContField2D::v_ImposeDirichletConditions for more detail.
677  map<int, vector<MultiRegions::ExtraDirDof> > &extraDirDofs =
678  m_assemblyMap->GetExtraDirDofs();
679  map<int, vector<MultiRegions::ExtraDirDof> >::iterator it;
680 
681  for (nv = 0; nv < nVel; ++nv)
682  {
684  = m_fields[nv]->GetBndCondExpansions();
685 
686  // First try to do parallel stuff
687  for (it = extraDirDofs.begin(); it != extraDirDofs.end(); ++it)
688  {
689  for (i = 0; i < it->second.size(); ++i)
690  {
691  inout[it->second.at(i).get<1>()*nVel + nv] =
692  bndCondExp[it->first]->GetCoeffs()[
693  it->second.at(i).get<0>()]*it->second.at(i).get<2>();
694  }
695  }
696  }
697 
698  m_assemblyMap->UniversalAssemble(inout);
699 
700  // Counter for the local Dirichlet boundary to global ordering.
701  int bndcnt = 0;
702 
703  // Now assemble local boundary contributions.
704  for (nv = 0; nv < nVel; ++nv)
705  {
707  = m_fields[nv]->GetBndCondExpansions();
708  const Array<OneD, const int> &bndMap
709  = m_assemblyMap->GetBndCondCoeffsToGlobalCoeffsMap();
710 
711  for (i = 0; i < bndCondExp.num_elements(); ++i)
712  {
713  if (m_fields[nv]->GetBndConditions()[i]->GetBoundaryConditionType()
715  {
716  const Array<OneD,const NekDouble> &bndCoeffs =
717  bndCondExp[i]->GetCoeffs();
718 
719  for (j = 0; j < bndCondExp[i]->GetNcoeffs(); ++j)
720  {
721  NekDouble sign =
722  m_assemblyMap->GetBndCondCoeffsToGlobalCoeffsSign(
723  bndcnt);
724  inout[bndMap[bndcnt++]] = sign * bndCoeffs[j];
725  }
726  }
727  else
728  {
729  bndcnt += bndCondExp[i]->GetNcoeffs();
730  }
731  }
732  }
733 
734  //
735  // -- Perform solve
736  //
737 
738  // Assemble forcing into the RHS.
739  m_assemblyMap->Assemble(forCoeffs, rhs);
740 
741  // Negate RHS to be consistent with matrix definition.
742  Vmath::Neg(rhs.num_elements(), rhs, 1);
743 
744  // Solve.
745  linSys->Solve(rhs, inout, m_assemblyMap);
746 
747  // Scatter the global ordering back to the alternate local ordering.
748  Array<OneD, NekDouble> tmp(nVel * nCoeffs);
749  m_assemblyMap->GlobalToLocal(inout, tmp);
750 
751  //
752  // -- Postprocess
753  //
754 
755  // Scatter back to field degrees of freedom
756  for (nv = 0; nv < nVel; ++nv)
757  {
758  for (i = 0; i < m_fields[nv]->GetExpSize(); ++i)
759  {
760  int nBnd = m_bmap[i].num_elements();
761  int nInt = m_imap[i].num_elements();
762  int offset = m_fields[nv]->GetCoeff_Offset(i);
763 
764  for (j = 0; j < nBnd; ++j)
765  {
766  m_fields[nv]->UpdateCoeffs()[offset+m_bmap[i][j]] =
767  tmp[nVel*offset + nv*nBnd + j];
768  }
769  for (j = 0; j < nInt; ++j)
770  {
771  m_fields[nv]->UpdateCoeffs()[offset+m_imap[i][j]] =
772  tmp[nVel*(offset + nBnd) + nv*nInt + j];
773  }
774  }
775  m_fields[nv]->BwdTrans(m_fields[nv]->GetCoeffs(),
776  m_fields[nv]->UpdatePhys());
777  }
778 }
779 
780 /**
781  * @brief Given a block matrix for an element, construct its static condensation
782  * matrices.
783  *
784  * This routine essentially duplicates the logic present in many of the
785  * LocalRegions matrix generation routines to construct the statically condensed
786  * equivalent of mat to pass to the GlobalLinSys solver.
787  *
788  * @param n Element/block number.
789  * @param exp Pointer to expansion.
790  * @param mat Block matrix containing matrix operator.
791  */
793  const int n,
796 {
797  int i, j, k, l;
798  const int nVel = mat.GetRows();
799  const int nB = exp->NumBndryCoeffs();
800  const int nI = exp->GetNcoeffs() - nB;
801  const int nBnd = exp->NumBndryCoeffs() * nVel;
802  const int nInt = exp->GetNcoeffs() * nVel - nBnd;
803  const MatrixStorage s = eFULL; // Maybe look into doing symmetric
804  // version of this?
805 
806  DNekMatSharedPtr A =
807  MemoryManager<DNekMat>::AllocateSharedPtr(nBnd, nBnd, 0.0, s);
808  DNekMatSharedPtr B =
809  MemoryManager<DNekMat>::AllocateSharedPtr(nBnd, nInt, 0.0, s);
810  DNekMatSharedPtr C =
811  MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nBnd, 0.0, s);
812  DNekMatSharedPtr D =
813  MemoryManager<DNekMat>::AllocateSharedPtr(nInt, nInt, 0.0, s);
814 
815  for (i = 0; i < nVel; ++i)
816  {
817  for (j = 0; j < nVel; ++j)
818  {
819  // Boundary-boundary and boundary-interior
820  for (k = 0; k < nB; ++k)
821  {
822  for (l = 0; l < nB; ++l)
823  {
824  (*A)(k + i*nB, l + j*nB) =
825  (*mat[i][j])(m_bmap[n][k], m_bmap[n][l]);
826  }
827 
828  for (l = 0; l < nI; ++l)
829  {
830  (*B)(k + i*nB, l + j*nI) =
831  (*mat[i][j])(m_bmap[n][k], m_imap[n][l]);
832  }
833  }
834 
835  // Interior-boundary / interior-interior
836  for (k = 0; k < nI; ++k)
837  {
838  for (l = 0; l < nB; ++l)
839  {
840  (*C)(k + i*nI, l + j*nB) =
841  (*mat[i][j])(m_imap[n][k], m_bmap[n][l]);
842  }
843 
844  for (l = 0; l < nI; ++l)
845  {
846  (*D)(k + i*nI, l + j*nI) =
847  (*mat[i][j])(m_imap[n][k], m_imap[n][l]);
848  }
849  }
850  }
851  }
852 
853  // Construct static condensation matrices.
854  D->Invert();
855  (*B) = (*B)*(*D);
856  (*A) = (*A) - (*B)*(*C);
857 
858  DNekScalMatSharedPtr tmp_mat;
859  m_schurCompl->SetBlock(
861  1.0, A));
862  m_BinvD ->SetBlock(
864  1.0, B));
865  m_C ->SetBlock(
867  1.0, C));
868  m_Dinv ->SetBlock(
870  1.0, D));
871 }
872 
873 /**
874  * @brief Construct a LaplacianIJ matrix for a given expansion.
875  *
876  * This routine constructs matrices whose entries contain the evaluation of
877  *
878  * \f[ \partial_{\tt k1} \phi_i \partial_{\tt k2} \phi_j \,dx \f]
879  */
881  const int k1,
882  const int k2,
883  const NekDouble scale,
885 {
886  const int nCoeffs = exp->GetNcoeffs();
887  const int nPhys = exp->GetTotPoints();
888  int i;
889 
891  nCoeffs, nCoeffs, 0.0, eFULL);
892 
893  Array<OneD, NekDouble> tmp2(nPhys);
894  Array<OneD, NekDouble> tmp3(nPhys);
895 
896  for (i = 0; i < nCoeffs; ++i)
897  {
898  Array<OneD, NekDouble> tmp1(nCoeffs, 0.0);
899  tmp1[i] = 1.0;
900 
901  exp->BwdTrans ( tmp1, tmp2);
902  exp->PhysDeriv (k1, tmp2, tmp3);
903  exp->IProductWRTDerivBase(k2, tmp3, tmp1);
904 
905  Vmath::Smul(
906  nCoeffs, scale, &tmp1[0], 1, &(ret->GetPtr())[0]+i*nCoeffs, 1);
907  }
908 
909  return ret;
910 }
911 
913  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
914  std::vector<std::string> &variables)
915 {
916  const int nVel = m_fields[0]->GetCoordim(0);
917  const int nCoeffs = m_fields[0]->GetNcoeffs();
918 
919  if (m_temperature.num_elements() == 0)
920  {
921  return;
922  }
923 
924  for (int i = 0; i < nVel; ++i)
925  {
926  Array<OneD, NekDouble> tFwd(nCoeffs);
927  m_fields[i]->FwdTrans(m_temperature[i], tFwd);
928  fieldcoeffs.push_back(tFwd);
929  variables.push_back(
930  "ThermStressDiv" + boost::lexical_cast<std::string>(i));
931  }
932 
933  if (m_stress.num_elements() == 0)
934  {
935  return;
936  }
937 
938  for (int i = 0; i < nVel; ++i)
939  {
940  for (int j = 0; j < nVel; ++j)
941  {
942  Array<OneD, NekDouble> tFwd(nCoeffs);
943  m_fields[i]->FwdTrans(m_stress[i][j], tFwd);
944  fieldcoeffs.push_back(tFwd);
945  variables.push_back(
946  "ThermStress"
947  + boost::lexical_cast<std::string>(i)
948  + boost::lexical_cast<std::string>(j));
949  }
950  }
951 }
952 
953 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_InitObject()
Set up the linear elasticity system.
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.
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
A base class for describing how to solve specific equations.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
NekDouble m_nu
Poisson ratio.
STL namespace.
virtual void v_DoSolve()
Solve elliptic linear elastic system.
NekDouble m_beta
Parameter dictating amount of thermal stress to add.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:293
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:251
DNekScalBlkMatSharedPtr m_BinvD
Matrix of elemental components.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, Array< OneD, NekDouble > > m_temperature
Storage for the temperature terms.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_stress
Storage for the thermal stress terms.
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Initialisation object for EquationSystem.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
static PreconditionerSharedPtr NullPreconditionerSharedPtr
This class is the abstraction of a global continuous two- dimensional spectral/hp element expansion w...
Definition: ContField2D.h:56
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Generate summary at runtime.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:226
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
DNekScalBlkMatSharedPtr m_Dinv
Matrix of elemental components.
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, Array< OneD, unsigned int > > m_imap
Interior maps for each of the fields.
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
Describe a linear system.
DNekScalBlkMatSharedPtr m_schurCompl
Schur complement boundary-boundary matrix.
EquationSystemFactory & GetEquationSystemFactory()
DNekScalBlkMatSharedPtr m_C
Matrix of elemental components.
NekDouble m_E
Young's modulus.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
DNekMatSharedPtr BuildLaplacianIJMatrix(const int k1, const int k2, const NekDouble scale, LocalRegions::ExpansionSharedPtr exp)
Construct a LaplacianIJ matrix for a given expansion.
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
CoupledAssemblyMapSharedPtr m_assemblyMap
Assembly map for the coupled (u,v,w) system.
Array< OneD, Array< OneD, unsigned int > > m_bmap
Boundary maps for each of the fields.
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:52
SOLVER_UTILS_EXPORT void SessionSummary(SummaryList &vSummary)
Write out a session summary.
void BuildMatrixSystem()
Build matrix system for linear elasticity equations.
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:191
DNekMat MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Geometry is curved or has non-constant factors.
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
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.cpp:285