MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_CoupledRBMFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2013 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_COUPLEDRBMFACTORY_DEF_HPP
47#define MUELU_COUPLEDRBMFACTORY_DEF_HPP
48
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_MultiVectorFactory.hpp>
51#include <Xpetra_VectorFactory.hpp>
52
54#include "MueLu_Level.hpp"
55#include "MueLu_Monitor.hpp"
56
57namespace MueLu {
58
59 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
61
62 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
64 if (currentLevel.IsAvailable(nspName_, NoFactory::get()) == false && currentLevel.GetLevelID() == 0) {
65 Input(currentLevel, "A");
66 //Input(currentLevel,"Coordinates");
67 }
68 if (currentLevel.GetLevelID() !=0) {
69 currentLevel.DeclareInput("Nullspace", GetFactory(nspName_).get(), this); /* ! "Nullspace" and nspName_ mismatch possible here */
70 }
71 }
72
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75 FactoryMonitor m(*this, "Structural acoustics nullspace factory", currentLevel);
76 RCP<MultiVector> nullspace;
77 if (currentLevel.GetLevelID() == 0) {
78 if (currentLevel.IsAvailable(nspName_, NoFactory::get())) {
79 nullspace = currentLevel.Get< RCP<MultiVector> >(nspName_, NoFactory::get());
80 GetOStream(Runtime1) << "Use user-given rigid body modes " << nspName_ << ": nullspace dimension=" << nullspace->getNumVectors() << " nullspace length=" << nullspace->getGlobalLength() << std::endl;
81 }
82 else {
83 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
84 RCP<MultiVector> Coords = Get< RCP<MultiVector> >(currentLevel,"Coordinates");
85 GetOStream(Runtime1) << "Generating nullspace for structural acoustics: dimension = " << numPDEs_ << std::endl;
86 RCP<const Map> xmap=A->getDomainMap();
87 nullspace = MultiVectorFactory::Build(xmap, 6);
88 Scalar zero = (Scalar) 0.0;
89 nullspace -> putScalar(zero);
90 ArrayRCP<Scalar> xnodes, ynodes, znodes;
91 Scalar cx, cy, cz;
92 ArrayRCP<Scalar> nsValues0, nsValues1, nsValues2, nsValues3, nsValues4, nsValues5;
93 int nDOFs=xmap->getLocalNumElements();
94 xnodes = Coords->getDataNonConst(0);
95 ynodes = Coords->getDataNonConst(1);
96 znodes = Coords->getDataNonConst(2);
97 cx = Coords->getVector(0)->meanValue();
98 cy = Coords->getVector(1)->meanValue();
99 cz = Coords->getVector(2)->meanValue();
100 nsValues0 = nullspace->getDataNonConst(0);
101 nsValues1 = nullspace->getDataNonConst(1);
102 nsValues2 = nullspace->getDataNonConst(2);
103 nsValues3 = nullspace->getDataNonConst(3);
104 nsValues4 = nullspace->getDataNonConst(4);
105 nsValues5 = nullspace->getDataNonConst(5);
106 for (int j=0; j<nDOFs; j+=numPDEs_) {
107 Scalar one = (Scalar) 1.0;
108 if( xmap->getGlobalElement(j) >= lastAcousticDOF_ ) {
109 Scalar xdiff = xnodes[j]-cx;
110 Scalar ydiff = ynodes[j]-cy;
111 Scalar zdiff = znodes[j]-cz;
112 // translation
113 nsValues0[j+0] = one;
114 nsValues1[j+1] = one;
115 nsValues2[j+2] = one;
116 // rotate around z-axis (x-y plane)
117 nsValues3[j+0] = -ydiff;
118 nsValues3[j+1] = xdiff;
119 // rotate around x-axis (y-z plane)
120 nsValues4[j+1] = -zdiff;
121 nsValues4[j+2] = ydiff;
122 // rotate around y-axis (x-z plane)
123 nsValues5[j+0] = zdiff;
124 nsValues5[j+2] = -xdiff;
125 }
126 else {
127 // translation
128 nsValues0[j+0] = one;
129 // insert random values and keep the top row for this node empty
130 nsValues1[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
131 nsValues1[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
132 nsValues2[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
133 nsValues2[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
134 nsValues3[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
135 nsValues3[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
136 nsValues4[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
137 nsValues4[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
138 nsValues5[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
139 nsValues5[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
140 }
141 }
142 } // end if "Nullspace" not available
143 }
144 else {
145 nullspace = currentLevel.Get< RCP<MultiVector> >("Nullspace", GetFactory(nspName_).get());
146 }
147 Set(currentLevel, "Nullspace", nullspace);
148 }
149
150 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
151 void CoupledRBMFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildRBM(RCP<Matrix>& A, RCP<MultiVector>& Coords, RCP<MultiVector>& nullspace) const {
152 GetOStream(Runtime1) << "Generating nullspace for structural acoustics: dimension = " << numPDEs_ << std::endl;
153 RCP<const Map> xmap=A->getDomainMap();
154 nullspace = MultiVectorFactory::Build(xmap, 6);
155 Scalar zero = (Scalar) 0.0;
156 nullspace -> putScalar(zero);
157 ArrayRCP<Scalar> xnodes, ynodes, znodes;
158 Scalar cx, cy, cz;
159 ArrayRCP<Scalar> nsValues0, nsValues1, nsValues2, nsValues3, nsValues4, nsValues5;
160 int nDOFs=xmap->getLocalNumElements();
161 xnodes = Coords->getDataNonConst(0);
162 ynodes = Coords->getDataNonConst(1);
163 znodes = Coords->getDataNonConst(2);
164 cx = Coords->getVector(0)->meanValue();
165 cy = Coords->getVector(1)->meanValue();
166 cz = Coords->getVector(2)->meanValue();
167 nsValues0 = nullspace->getDataNonConst(0);
168 nsValues1 = nullspace->getDataNonConst(1);
169 nsValues2 = nullspace->getDataNonConst(2);
170 nsValues3 = nullspace->getDataNonConst(3);
171 nsValues4 = nullspace->getDataNonConst(4);
172 nsValues5 = nullspace->getDataNonConst(5);
173 for (int j=0; j<nDOFs; j+=numPDEs_) {
174 Scalar one = (Scalar) 1.0;
175 if( xmap->getGlobalElement(j) >= lastAcousticDOF_ ) {
176 Scalar xdiff = xnodes[j]-cx;
177 Scalar ydiff = ynodes[j]-cy;
178 Scalar zdiff = znodes[j]-cz;
179 // translation
180 nsValues0[j+0] = one;
181 nsValues1[j+1] = one;
182 nsValues2[j+2] = one;
183 // rotate around z-axis (x-y plane)
184 nsValues3[j+0] = -ydiff;
185 nsValues3[j+1] = xdiff;
186 // rotate around x-axis (y-z plane)
187 nsValues4[j+1] = -zdiff;
188 nsValues4[j+2] = ydiff;
189 // rotate around y-axis (x-z plane)
190 nsValues5[j+0] = zdiff;
191 nsValues5[j+2] = -xdiff;
192 }
193 else {
194 // translation
195 nsValues0[j+0] = one;
196 // insert random values and keep the top row for this node empty
197 nsValues1[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
198 nsValues1[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
199 nsValues2[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
200 nsValues2[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
201 nsValues3[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
202 nsValues3[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
203 nsValues4[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
204 nsValues4[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
205 nsValues5[j+1] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
206 nsValues5[j+2] = (Scalar) (((double) rand()) / ((double) RAND_MAX));
207 }
208 }
209 }
210
211} // namespace MueLu
212
213#define MUELU_COUPLEDRBMFACTORY_SHORT
214#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.
Definition: MueLu_Level.hpp:99
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.
Definition: MueLu_Level.cpp:76
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)