Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CoupledAssemblyMap.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: CoupledAssemblyMap.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: Coupled assembly map for linear elasticity solver.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 #include <LocalRegions/SegExp.h>
42 
43 #include <iomanip>
44 
45 using namespace std;
46 
47 namespace Nektar
48 {
49 
51 
52 /**
53  * @brief Take an existing assembly map and create a coupled version suitable
54  * for use in the linear elasticity solver.
55  *
56  * The linear elasticity solver requires a slight reordering of local and global
57  * coefficients to support problems of the form
58  *
59  * [ A B C ] [ u ] = [ f_u ]
60  * [ D E F ] [ v ] [ f_v ]
61  * [ G H I ] [ w ] [ f_w ]
62  */
63 
64 CoupledAssemblyMap::CoupledAssemblyMap(
68  const Array<OneD, const BoundaryCondShPtr> &boundaryConditions,
70  AssemblyMapCG(pSession)
71 {
72  int nVel = fields[0]->GetCoordim(0);
73 
74  // Multi-level static condensation doesn't work yet.
78  "Multi-level static condensation not supported.");
79 
80  // Copy various coefficient counts, and multiply by the dimension of the
81  // problem to obtain our new values.
82  m_numLocalDirBndCoeffs = cgMap->GetNumLocalDirBndCoeffs() * nVel;
83  m_numLocalBndCoeffs = cgMap->GetNumLocalBndCoeffs() * nVel;
84  m_numLocalCoeffs = cgMap->GetNumLocalCoeffs() * nVel;
85  m_numGlobalBndCoeffs = cgMap->GetNumGlobalBndCoeffs() * nVel;
86  m_numGlobalDirBndCoeffs = cgMap->GetNumGlobalDirBndCoeffs() * nVel;
87  m_numGlobalCoeffs = cgMap->GetNumGlobalCoeffs() * nVel;
88  m_signChange = cgMap->GetSignChange();
89  m_systemSingular = cgMap->GetSingularSystem();
90 
91  // Copy static condensation information. TODO: boundary and interior patches
92  // need to be re-ordered in order to allow for multi-level static
93  // condensation support.
94  m_staticCondLevel = cgMap->GetStaticCondLevel();
95  m_lowestStaticCondLevel = cgMap->GetLowestStaticCondLevel();
96  m_numPatches = cgMap->GetNumPatches();
97  m_numLocalBndCoeffsPerPatch = cgMap->GetNumLocalBndCoeffsPerPatch();
98  m_numLocalIntCoeffsPerPatch = cgMap->GetNumLocalIntCoeffsPerPatch();
99 
100  // Set up local to global and boundary condition maps.
101  const int nLocBndCondDofs = cgMap->
102  GetBndCondCoeffsToGlobalCoeffsMap().num_elements() * nVel;
103 
104  // Allocate storage for local to global maps.
110  Array<OneD, int>(nLocBndCondDofs,-1);
111 
112  // Only require a sign map if we are using modal polynomials in the
113  // expansion and the order is >= 3.
114  if (m_signChange)
115  {
121  Array<OneD, NekDouble>(nLocBndCondDofs,1.0);
122  }
123  else
124  {
128  }
129 
130  const LocalRegions::ExpansionVector &locExpVector
131  = *(fields[0]->GetExp());
132 
133  map<int, int> newGlobalIds;
134  int i, j, n, cnt1, cnt2;
135 
136  // Order local boundary degrees of freedom. These are basically fine; we
137  // reorder storage so that we loop over each element and then each component
138  // of velocity, by applying a mapping l2g -> nVel*l2g + n, for 0 <= n <
139  // nVel. Note that Dirichlet ordering is preserved under this
140  // transformation.
141  cnt1 = cnt2 = 0;
142  for (i = 0; i < locExpVector.size(); ++i)
143  {
144  const int nBndCoeffs = locExpVector[i]->NumBndryCoeffs();
145 
146  for (n = 0; n < nVel; ++n)
147  {
148  for (j = 0; j < nBndCoeffs; ++j, ++cnt1)
149  {
150  const int l2g = cgMap->GetLocalToGlobalBndMap()[cnt2+j];
151  m_localToGlobalBndMap[cnt1] = nVel * l2g + n;
152 
153  if (m_signChange)
154  {
155  m_localToGlobalBndSign[cnt1] =
156  cgMap->GetLocalToGlobalBndSign()[cnt2+j];
157  }
158 
159  if (n == 0)
160  {
161  const int l2gnew = m_localToGlobalBndMap[cnt1];
162  if (newGlobalIds.count(l2g))
163  {
164  ASSERTL1(newGlobalIds[l2g] == l2gnew,
165  "Consistency error");
166  }
167  newGlobalIds[l2g] = l2gnew;
168  }
169  }
170  }
171 
172  cnt2 += nBndCoeffs;
173  }
174 
175  // Grab map of extra Dirichlet degrees of freedom for parallel runs
176  // (particularly in 3D).
177  m_extraDirDofs = cgMap->GetExtraDirDofs();
178 
179  // Counter for remaining interior degrees of freedom.
180  int globalId = m_numGlobalBndCoeffs;
181 
182  // Interior degrees of freedom are a bit more tricky -- global linear system
183  // solve relies on them being in the same order as the BinvD, C and invD
184  // matrices.
185  cnt1 = cnt2 = 0;
186  for (i = 0; i < locExpVector.size(); ++i)
187  {
188  const int nCoeffs = locExpVector[i]->GetNcoeffs();
189  const int nBndCoeffs = locExpVector[i]->NumBndryCoeffs();
190 
191  for (n = 0; n < nVel; ++n)
192  {
193  for (j = 0; j < nBndCoeffs; ++j, ++cnt1, ++cnt2)
194  {
195  const int l2g = m_localToGlobalBndMap[cnt2];
196  m_localToGlobalMap[cnt1] = l2g;
197 
198  if (m_signChange)
199  {
201  }
202  }
203  }
204 
205  for (n = 0; n < nVel; ++n)
206  {
207  for (j = 0; j < nCoeffs - nBndCoeffs; ++j, ++cnt1)
208  {
209  m_localToGlobalMap[cnt1] = globalId++;
210  }
211  }
212  }
213 
214  for (i = 0; i < m_localToGlobalMap.num_elements(); ++i)
215  {
216  ASSERTL1(m_localToGlobalMap[i] != -1, "Consistency error");
217  }
218 
219  ASSERTL1(globalId == m_numGlobalCoeffs, "Consistency error");
220 
221  cnt1 = 0;
222  for (n = 0; n < nVel; ++n)
223  {
224  for (i = 0; i < nLocBndCondDofs/nVel; ++i, ++cnt1)
225  {
226  const int l2g = cgMap->GetBndCondCoeffsToGlobalCoeffsMap()[i];
227  int newId = newGlobalIds[l2g];
228  m_bndCondCoeffsToGlobalCoeffsMap[cnt1] = newId + n;
229 
230  if (m_signChange)
231  {
233  cgMap->GetBndCondCoeffsToGlobalCoeffsSign(i);
234  }
235  }
236  }
237 
238  // Finally, set up global to universal maps.
243 
244  for (i = 0; i < cgMap->GetNumGlobalBndCoeffs(); ++i)
245  {
246  for (n = 0; n < nVel; ++n)
247  {
248  m_globalToUniversalBndMap[i*nVel + n] =
249  cgMap->GetGlobalToUniversalBndMap()[i]*nVel + n;
250  m_globalToUniversalMap[i*nVel + n] =
251  cgMap->GetGlobalToUniversalBndMap()[i]*nVel + n;
252  }
253  }
254 
258  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
259  {
260  tmp[i] = m_globalToUniversalBndMap[i];
261  }
262 
263  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
264  m_gsh = Gs::Init(tmp, vCommRow);
265  m_bndGsh = Gs::Init(tmp2, vCommRow);
266  Gs::Unique(tmp, vCommRow);
267  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
268  {
269  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
270  }
271  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
272  {
273  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
274  }
275 
276  m_hash = boost::hash_range(
277  m_localToGlobalMap.begin(), m_localToGlobalMap.end());
278 }
279 
280 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:320
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:314
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:306
static Array< OneD, NekDouble > NullNekDouble1DArray
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
STL namespace.
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:331
std::map< int, std::vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:309
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:388
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:363
Array< OneD, NekDouble > m_bndCondCoeffsToGlobalCoeffsSign
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:354
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
static void Unique(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Updates pId to negate all-but-one references to each universal ID.
Definition: GsLib.hpp:184
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:390
int m_lowestStaticCondLevel
Lowest static condensation level.
Definition: AssemblyMap.h:397
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
Definition: AssemblyMap.h:348
const Array< OneD, const int > & GetBndCondCoeffsToGlobalCoeffsMap()
Retrieves the global indices corresponding to the boundary expansion modes.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:316
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:352
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:312
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:384
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:350
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:358
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:360
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:208
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:52
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:386