Nektar++
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// 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: Coupled assembly map for linear elasticity solver.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36
41#include <LocalRegions/SegExp.h>
44
45namespace Nektar
46{
47
49
50/**
51 * @brief Take an existing assembly map and create a coupled version suitable
52 * for use in the linear elasticity solver.
53 *
54 * The linear elasticity solver requires a slight reordering of local and global
55 * coefficients to support problems of the form
56 *
57 * [ A B C ] [ u ] = [ f_u ]
58 * [ D E F ] [ v ] [ f_v ]
59 * [ G H I ] [ w ] [ f_w ]
60 */
61
64 [[maybe_unused]] const SpatialDomains::MeshGraphSharedPtr &graph,
66 [[maybe_unused]] const Array<OneD, const BoundaryCondShPtr>
67 &boundaryConditions,
69 : AssemblyMapCG(pSession, fields[0]->GetComm())
70{
71 int nVel = fields[0]->GetCoordim(0);
72
73 // Multi-level static condensation doesn't work yet.
77 "Multi-level static condensation not supported.");
78
79 // Copy various coefficient counts, and multiply by the dimension of the
80 // problem to obtain our new values.
81 m_numLocalDirBndCoeffs = cgMap->GetNumLocalDirBndCoeffs() * nVel;
82 m_numLocalBndCoeffs = cgMap->GetNumLocalBndCoeffs() * nVel;
83 m_numLocalCoeffs = cgMap->GetNumLocalCoeffs() * nVel;
84 m_numGlobalBndCoeffs = cgMap->GetNumGlobalBndCoeffs() * nVel;
85 m_numGlobalDirBndCoeffs = cgMap->GetNumGlobalDirBndCoeffs() * nVel;
86 m_numGlobalCoeffs = cgMap->GetNumGlobalCoeffs() * nVel;
87 m_signChange = cgMap->GetSignChange();
88 m_systemSingular = cgMap->GetSingularSystem();
89
90 // Copy static condensation information. TODO: boundary and interior patches
91 // need to be re-ordered in order to allow for multi-level static
92 // condensation support.
93 m_staticCondLevel = cgMap->GetStaticCondLevel();
94 m_lowestStaticCondLevel = cgMap->GetLowestStaticCondLevel();
95 m_numPatches = cgMap->GetNumPatches();
96 m_numLocalBndCoeffsPerPatch = cgMap->GetNumLocalBndCoeffsPerPatch();
97 m_numLocalIntCoeffsPerPatch = cgMap->GetNumLocalIntCoeffsPerPatch();
98
99 // Allocate storage for local to global maps.
105
106 // Only require a sign map if we are using modal polynomials in the
107 // expansion and the order is >= 3.
108 if (m_signChange)
109 {
113 }
114 else
115 {
118 }
119
120 const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
121
122 std::map<int, int> newGlobalIds;
123 int i, j, n, cnt1, cnt2;
124
125 // Order local boundary degrees of freedom. These are basically fine; we
126 // reorder storage so that we loop over each element and then each component
127 // of velocity, by applying a mapping l2g -> nVel*l2g + n, for 0 <= n <
128 // nVel. Note that Dirichlet ordering is preserved under this
129 // transformation.
130 cnt1 = cnt2 = 0;
131 for (i = 0; i < locExpVector.size(); ++i)
132 {
133 const int nBndCoeffs = locExpVector[i]->NumBndryCoeffs();
134
135 for (n = 0; n < nVel; ++n)
136 {
137
138 for (j = 0; j < nBndCoeffs; ++j, ++cnt1)
139 {
140 const int l2g = cgMap->GetLocalToGlobalBndMap()[cnt2 + j];
141 m_localToGlobalBndMap[cnt1] = nVel * l2g + n;
142
143 if (m_signChange)
144 {
146 cgMap->GetLocalToGlobalBndSign()[cnt2 + j];
147 }
148
149 if (n == 0)
150 {
151 const int l2gnew = m_localToGlobalBndMap[cnt1];
152 if (newGlobalIds.count(l2g))
153 {
154 ASSERTL1(newGlobalIds[l2g] == l2gnew,
155 "Consistency error");
156 }
157 newGlobalIds[l2g] = l2gnew;
158 }
159 }
160 }
161
162 cnt2 += nBndCoeffs;
163 }
164
165 // Counter for remaining interior degrees of freedom.
166 int globalId = m_numGlobalBndCoeffs;
167
168 // Interior degrees of freedom are a bit more tricky -- global linear system
169 // solve relies on them being in the same order as the BinvD, C and invD
170 // matrices.
171
172 // Also set up the localToBndMap and localTolocalIntMap which just
173 // take out the boundary blocks and interior blocks from the input
174 // ordering where we have bnd and interior for each elements
175 cnt1 = cnt2 = 0;
176 int bnd_cnt = 0;
177 int int_cnt = 0;
178 for (i = 0; i < locExpVector.size(); ++i)
179 {
180 const int nCoeffs = locExpVector[i]->GetNcoeffs();
181 const int nBndCoeffs = locExpVector[i]->NumBndryCoeffs();
182
183 for (n = 0; n < nVel; ++n)
184 {
185 for (j = 0; j < nBndCoeffs; ++j, ++cnt1, ++cnt2)
186 {
187 m_localToLocalBndMap[bnd_cnt++] = cnt1;
188
189 const int l2g = m_localToGlobalBndMap[cnt2];
190 m_localToGlobalMap[cnt1] = l2g;
191
192 if (m_signChange)
193 {
195 }
196 }
197 }
198
199 for (n = 0; n < nVel; ++n)
200 {
201 for (j = 0; j < nCoeffs - nBndCoeffs; ++j, ++cnt1)
202 {
203 m_localToLocalIntMap[int_cnt++] = cnt1;
204
205 m_localToGlobalMap[cnt1] = globalId++;
206 }
207 }
208 }
209
210 for (i = 0; i < m_localToGlobalMap.size(); ++i)
211 {
212 ASSERTL1(m_localToGlobalMap[i] != -1, "Consistency error");
213 }
214
215 ASSERTL1(globalId == m_numGlobalCoeffs, "Consistency error");
216
217 // Finally, set up global to universal maps.
222
223 for (i = 0; i < cgMap->GetNumGlobalBndCoeffs(); ++i)
224 {
225 for (n = 0; n < nVel; ++n)
226 {
227 m_globalToUniversalBndMap[i * nVel + n] =
228 cgMap->GetGlobalToUniversalBndMap()[i] * nVel + n;
229 m_globalToUniversalMap[i * nVel + n] =
230 cgMap->GetGlobalToUniversalBndMap()[i] * nVel + n;
231 }
232 }
233
237 for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
238 {
239 tmp[i] = m_globalToUniversalBndMap[i];
240 }
241
242 LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
243 m_gsh = Gs::Init(tmp, vCommRow);
244 m_bndGsh = Gs::Init(tmp2, vCommRow);
245 Gs::Unique(tmp, vCommRow);
246 for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
247 {
248 m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
249 }
250 for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
251 {
252 m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
253 }
254
256}
257
258} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
CoupledAssemblyMap(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const MultiRegions::AssemblyMapCGSharedPtr &cgMap, const Array< OneD, const BoundaryCondShPtr > &boundaryConditions, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Take an existing assembly map and create a coupled version suitable for use in the linear elasticity ...
Constructs mappings for the C0 scalar continuous Galerkin formulation.
Definition: AssemblyMapCG.h:66
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
int m_lowestStaticCondLevel
Lowest static condensation level.
Definition: AssemblyMap.h:436
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:399
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:361
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:375
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
Definition: AssemblyMap.h:384
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:372
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:394
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:380
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:427
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:342
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:423
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:346
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:348
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:429
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:350
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
Definition: AssemblyMap.h:378
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:396
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:333
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
Definition: AssemblyMap.h:382
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:425
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:344
static gs_data * Init(const Nektar::Array< OneD, long > &pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:190
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:225
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:68
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:50
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:213
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::size_t hash_range(Iter first, Iter last)
Definition: HashUtils.hpp:64
static Array< OneD, NekDouble > NullNekDouble1DArray
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273