MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_CoupledRBMFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MueLu: A package for multigrid based preconditioning
4//
5// Copyright 2012 NTESS and the MueLu contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef MUELU_COUPLEDRBMFACTORY_DEF_HPP
11#define MUELU_COUPLEDRBMFACTORY_DEF_HPP
12
13#include <Xpetra_Matrix.hpp>
14#include <Xpetra_MultiVectorFactory.hpp>
15#include <Xpetra_VectorFactory.hpp>
16
18#include "MueLu_Level.hpp"
19#include "MueLu_Monitor.hpp"
20
21namespace MueLu {
22
23template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
25
26template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
28 if (currentLevel.IsAvailable(nspName_, NoFactory::get()) == false && currentLevel.GetLevelID() == 0) {
29 Input(currentLevel, "A");
30 // Input(currentLevel,"Coordinates");
31 }
32 if (currentLevel.GetLevelID() != 0) {
33 currentLevel.DeclareInput("Nullspace", GetFactory(nspName_).get(), this); /* ! "Nullspace" and nspName_ mismatch possible here */
34 }
35}
36
37template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
39 FactoryMonitor m(*this, "Structural acoustics nullspace factory", currentLevel);
40 RCP<MultiVector> nullspace;
41 if (currentLevel.GetLevelID() == 0) {
42 if (currentLevel.IsAvailable(nspName_, NoFactory::get())) {
43 nullspace = currentLevel.Get<RCP<MultiVector> >(nspName_, NoFactory::get());
44 GetOStream(Runtime1) << "Use user-given rigid body modes " << nspName_ << ": nullspace dimension=" << nullspace->getNumVectors() << " nullspace length=" << nullspace->getGlobalLength() << std::endl;
45 } else {
46 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel, "A");
47 RCP<MultiVector> Coords = Get<RCP<MultiVector> >(currentLevel, "Coordinates");
48 GetOStream(Runtime1) << "Generating nullspace for structural acoustics: dimension = " << numPDEs_ << std::endl;
49 RCP<const Map> xmap = A->getDomainMap();
50 nullspace = MultiVectorFactory::Build(xmap, 6);
51 Scalar zero = (Scalar)0.0;
52 nullspace->putScalar(zero);
53 ArrayRCP<Scalar> xnodes, ynodes, znodes;
54 Scalar cx, cy, cz;
55 ArrayRCP<Scalar> nsValues0, nsValues1, nsValues2, nsValues3, nsValues4, nsValues5;
56 int nDOFs = xmap->getLocalNumElements();
57 xnodes = Coords->getDataNonConst(0);
58 ynodes = Coords->getDataNonConst(1);
59 znodes = Coords->getDataNonConst(2);
60 cx = Coords->getVector(0)->meanValue();
61 cy = Coords->getVector(1)->meanValue();
62 cz = Coords->getVector(2)->meanValue();
63 nsValues0 = nullspace->getDataNonConst(0);
64 nsValues1 = nullspace->getDataNonConst(1);
65 nsValues2 = nullspace->getDataNonConst(2);
66 nsValues3 = nullspace->getDataNonConst(3);
67 nsValues4 = nullspace->getDataNonConst(4);
68 nsValues5 = nullspace->getDataNonConst(5);
69 for (int j = 0; j < nDOFs; j += numPDEs_) {
70 Scalar one = (Scalar)1.0;
71 if (xmap->getGlobalElement(j) >= lastAcousticDOF_) {
72 Scalar xdiff = xnodes[j] - cx;
73 Scalar ydiff = ynodes[j] - cy;
74 Scalar zdiff = znodes[j] - cz;
75 // translation
76 nsValues0[j + 0] = one;
77 nsValues1[j + 1] = one;
78 nsValues2[j + 2] = one;
79 // rotate around z-axis (x-y plane)
80 nsValues3[j + 0] = -ydiff;
81 nsValues3[j + 1] = xdiff;
82 // rotate around x-axis (y-z plane)
83 nsValues4[j + 1] = -zdiff;
84 nsValues4[j + 2] = ydiff;
85 // rotate around y-axis (x-z plane)
86 nsValues5[j + 0] = zdiff;
87 nsValues5[j + 2] = -xdiff;
88 } else {
89 // translation
90 nsValues0[j + 0] = one;
91 // insert random values and keep the top row for this node empty
92 nsValues1[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
93 nsValues1[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
94 nsValues2[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
95 nsValues2[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
96 nsValues3[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
97 nsValues3[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
98 nsValues4[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
99 nsValues4[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
100 nsValues5[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
101 nsValues5[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
102 }
103 }
104 } // end if "Nullspace" not available
105 } else {
106 nullspace = currentLevel.Get<RCP<MultiVector> >("Nullspace", GetFactory(nspName_).get());
107 }
108 Set(currentLevel, "Nullspace", nullspace);
109}
110
111template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112void CoupledRBMFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildRBM(RCP<Matrix>& A, RCP<MultiVector>& Coords, RCP<MultiVector>& nullspace) const {
113 GetOStream(Runtime1) << "Generating nullspace for structural acoustics: dimension = " << numPDEs_ << std::endl;
114 RCP<const Map> xmap = A->getDomainMap();
115 nullspace = MultiVectorFactory::Build(xmap, 6);
116 Scalar zero = (Scalar)0.0;
117 nullspace->putScalar(zero);
118 ArrayRCP<Scalar> xnodes, ynodes, znodes;
119 Scalar cx, cy, cz;
120 ArrayRCP<Scalar> nsValues0, nsValues1, nsValues2, nsValues3, nsValues4, nsValues5;
121 int nDOFs = xmap->getLocalNumElements();
122 xnodes = Coords->getDataNonConst(0);
123 ynodes = Coords->getDataNonConst(1);
124 znodes = Coords->getDataNonConst(2);
125 cx = Coords->getVector(0)->meanValue();
126 cy = Coords->getVector(1)->meanValue();
127 cz = Coords->getVector(2)->meanValue();
128 nsValues0 = nullspace->getDataNonConst(0);
129 nsValues1 = nullspace->getDataNonConst(1);
130 nsValues2 = nullspace->getDataNonConst(2);
131 nsValues3 = nullspace->getDataNonConst(3);
132 nsValues4 = nullspace->getDataNonConst(4);
133 nsValues5 = nullspace->getDataNonConst(5);
134 for (int j = 0; j < nDOFs; j += numPDEs_) {
135 Scalar one = (Scalar)1.0;
136 if (xmap->getGlobalElement(j) >= lastAcousticDOF_) {
137 Scalar xdiff = xnodes[j] - cx;
138 Scalar ydiff = ynodes[j] - cy;
139 Scalar zdiff = znodes[j] - cz;
140 // translation
141 nsValues0[j + 0] = one;
142 nsValues1[j + 1] = one;
143 nsValues2[j + 2] = one;
144 // rotate around z-axis (x-y plane)
145 nsValues3[j + 0] = -ydiff;
146 nsValues3[j + 1] = xdiff;
147 // rotate around x-axis (y-z plane)
148 nsValues4[j + 1] = -zdiff;
149 nsValues4[j + 2] = ydiff;
150 // rotate around y-axis (x-z plane)
151 nsValues5[j + 0] = zdiff;
152 nsValues5[j + 2] = -xdiff;
153 } else {
154 // translation
155 nsValues0[j + 0] = one;
156 // insert random values and keep the top row for this node empty
157 nsValues1[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
158 nsValues1[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
159 nsValues2[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
160 nsValues2[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
161 nsValues3[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
162 nsValues3[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
163 nsValues4[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
164 nsValues4[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
165 nsValues5[j + 1] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
166 nsValues5[j + 2] = (Scalar)(((double)rand()) / ((double)RAND_MAX));
167 }
168 }
169}
170
171} // namespace MueLu
172
173#define MUELU_COUPLEDRBMFACTORY_SHORT
174#endif // MUELU_COUPLEDRBMFACTORY_DEF_HPP
MueLu::DefaultScalar Scalar
void BuildRBM(RCP< Matrix > &A, RCP< MultiVector > &Coords, RCP< MultiVector > &nullspace) const
void Build(Level &currentLevel) const
Build an object with this factory.
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose)