MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
53 #include "MueLu_CoupledRBMFactory_decl.hpp"
54 #include "MueLu_Level.hpp"
55 #include "MueLu_Monitor.hpp"
56 
57 namespace MueLu {
58 
59  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
60  CoupledRBMFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~CoupledRBMFactory() {}
61 
62  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63  void CoupledRBMFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &currentLevel) const {
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>
74  void CoupledRBMFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const {
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->getNodeNumElements();
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->getNodeNumElements();
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::Level::Get
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
Definition: MueLu_Level.hpp:191
MueLu::CoupledRBMFactory::DeclareInput
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Definition: MueLu_CoupledRBMFactory_def.hpp:63
MueLu::CoupledRBMFactory::Build
void Build(Level &currentLevel) const
Build an object with this factory.
Definition: MueLu_CoupledRBMFactory_def.hpp:74
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
MueLu_Level.hpp
Xpetra_Matrix.hpp
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
Xpetra_MultiVectorFactory.hpp
MueLu::CoupledRBMFactory::BuildRBM
void BuildRBM(RCP< Matrix > &A, RCP< MultiVector > &Coords, RCP< MultiVector > &nullspace) const
Definition: MueLu_CoupledRBMFactory_def.hpp:151
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::CoupledRBMFactory::~CoupledRBMFactory
virtual ~CoupledRBMFactory()
Destructor.
Definition: MueLu_CoupledRBMFactory_def.hpp:60
MueLu_Monitor.hpp
Xpetra_VectorFactory.hpp
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:63
MueLu::Level::DeclareInput
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
Definition: MueLu_Level.cpp:150
MueLu_CoupledRBMFactory_decl.hpp
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Definition: MueLu_Level.hpp:332