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
std::string nspName_
name of nullspace vector on finest level
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.
void Input(Level &level, const std::string &varName) const
T Get(Level &level, const std::string &varName) const
void Set(Level &level, const std::string &varName, const T &data) const
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory().
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()
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose).