FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
beam.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
44//
45// This is a simple program to exercise FEI classes for the
46// purposes of testing, code tuning and scaling studies.
47//
48
49#include "fei_iostream.hpp"
50
51//Including the header fei_base.hpp pulls in the declaration for
52//various classes in the 'fei' namespace.
53
54#include "fei_base.hpp"
55
56
57//Make provision for using any one of several solver libraries. This is
58//handled by the code in LibraryFactory.{hpp,cpp}.
59
61
62
63//we need to include some 'test_utils' headers which are used for
64//setting up the data for the hex-beam example problem.
65
68
69
70//fei_test_utils.hpp declares things like the intialize_mpi function.
71
73
74//
75//Include definitions of macros like 'CHK_ERR' to call functions and check
76//the return code.
77//
78#undef fei_file
79#define fei_file "beam.cpp"
80#include "fei_ErrMacros.hpp"
81
82//==============================================================================
83//Here's the main...
84//==============================================================================
85int main(int argc, char** argv)
86{
87 MPI_Comm comm;
88 int numProcs, localProc;
89 CHK_ERR( fei_test_utils::initialize_mpi(argc, argv, localProc, numProcs) );
90 comm = MPI_COMM_WORLD;
91
92 double start_time = fei::utils::cpu_time();
93
94 //read input parameters from a file specified on the command-line with
95 // '-i file'
96 std::vector<std::string> stdstrings;
98 comm, localProc,
99 stdstrings) );
100 fei::ParameterSet params;
101 //parse the strings from the input file into a fei::ParameterSet object.
102 fei::utils::parse_strings(stdstrings, " ", params);
103
104 std::string solverName;
105 std::string datasource;
106 std::string constraintform;
107 int W = 0;
108 int D = 0;
109 int DofPerNode = 0;
110
111 int errcode = 0;
112 errcode += params.getStringParamValue("SOLVER_LIBRARY", solverName);
113 errcode += params.getStringParamValue("DATA_SOURCE", datasource);
114 errcode += params.getIntParamValue("W", W);
115 errcode += params.getIntParamValue("D", D);
116 errcode += params.getIntParamValue("DofPerNode", DofPerNode);
117
118 if (errcode != 0) {
119 fei::console_out() << "Failed to find one or more required parameters in input-file."
120 << std::endl << "Required parameters:"<<std::endl
121 << "SOLVER_LIBRARY" << std::endl
122 << "DATA_SOURCE" << std::endl
123 << "CONSTRAINT_FORM" << std::endl
124 << "W" << std::endl << "D" << std::endl << "DofPerNode" << std::endl;
125
126#ifndef FEI_SER
127 MPI_Finalize();
128#endif
129
130 return(-1);
131 }
132
133 params.getStringParamValue("CONSTRAINT_FORM", constraintform);
134
135 bool slave_constraints = (constraintform == "slaves");
136
137 HexBeam* hexcubeptr = NULL;
138 if (datasource == "HexBeam") {
139 hexcubeptr = new HexBeam(W, D, DofPerNode,
140 HexBeam::OneD, numProcs, localProc);
141 }
142 else {
143 hexcubeptr = new HexBeamCR(W, D, DofPerNode,
144 HexBeam::OneD, numProcs, localProc);
145 }
146
147 HexBeam& hexcube = *hexcubeptr;
148
149 if (localProc == 0) {
150 int numCRs = (W+1)*(W+1)* ((numProcs*2)-1);
151 if (hexcube.getNumCRs() < 1) numCRs = 0;
152 //macros std::cout and std::endl are aliases for std::cout and std::endl,
153 //defined in fei_iostream.hpp.
154 std::cout << std::endl;
155 std::cout << "========================================================\n";
156 std::cout << "FEI version: " << fei::utils::version() << std::endl;
157 std::cout << "--------------------------------------------------------\n";
158 std::cout << "Size W: " << W << " (num-elements-along-side-of-cube)\n";
159 std::cout << "Size D: " << D << " (num-elements-along-depth-of-cube)\n";
160 std::cout << "DOF per node: " << DofPerNode <<"\n";
161 std::cout << "Num local elements: " << hexcube.localNumElems_ << "\n";
162 std::cout << "Num global elements: " << hexcube.totalNumElems_ << "\n";
163 std::cout << "Num local DOF: " << hexcube.numLocalDOF_ << "\n";
164 std::cout << "Num global DOF: " << hexcube.numGlobalDOF_ << "\n";
165 std::cout << "Num global CRs: " << numCRs << "\n";
166 std::cout << "========================================================"
167 << std::endl;
168 }
169
170 double start_init_time = fei::utils::cpu_time();
171
173 try {
174 factory = fei::create_fei_Factory(comm, solverName.c_str());
175 }
176 catch (...) {
177 std::cout << "library " << solverName << " not available."<<std::endl;
178
179#ifndef FEI_SER
180 MPI_Finalize();
181#endif
182
183 return(-1);
184 }
185 if (factory.get() == NULL) {
186 std::cout << "fei::Factory allocation failed." << std::endl;
187
188#ifndef FEI_SER
189 MPI_Finalize();
190#endif
191
192 return(-1);
193 }
194
195 factory->parameters(params);
196
198 factory->createVectorSpace(comm, "beam");
199
202 factory->createMatrixGraph(nodeSpace, dummy, "beam");
203
204 //load some solver-control parameters.
205 nodeSpace->setParameters(params);
206 matrixGraph->setParameters(params);
207
208 int fieldID = 0;
209 int fieldSize = hexcube.numDofPerNode();
210 int nodeIDType = 0;
211 int crIDType = 1;
212
213 nodeSpace->defineFields( 1, &fieldID, &fieldSize );
214 nodeSpace->defineIDTypes(1, &nodeIDType );
215 nodeSpace->defineIDTypes(1, &crIDType );
216
217 //HexBeam_Functions::init_elem_connectivities demonstrates the process
218 //of initializing element-to-node connectivities on the fei::MatrixGraph.
220 init_elem_connectivities( matrixGraph.get(), hexcube ) );
221
222 //HexBeam_Functions::init_shared_nodes demonstrates the process of
223 //initializing shared nodes on the fei::MatrixGraph.
225 init_shared_nodes( matrixGraph.get(), hexcube ) );
226
227 int firstLocalCRID = 0;
228 if (slave_constraints) {
230 init_slave_constraints( matrixGraph.get(), hexcube) );
231 }
232 else {
234 init_constraints( matrixGraph.get(), hexcube, localProc, firstLocalCRID) );
235 }
236
237 CHK_ERR( matrixGraph->initComplete() );
238
239 double fei_init_time = fei::utils::cpu_time() - start_init_time;
240
241 if (localProc == 0) {
242 std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
243 std::cout << "Initialization cpu time: " << fei_init_time << std::endl;
244 }
245
246 //Now the initialization phase is complete. Next we'll do the load phase,
247 //which for this problem just consists of loading the element data
248 //(element-wise stiffness arrays and load vectors) and the boundary
249 //condition data.
250
251 double fei_creatematrix_start_time = fei::utils::cpu_time();
252
253 fei::SharedPtr<fei::Matrix> mat = factory->createMatrix(matrixGraph);
254
255 double fei_creatematrix_time = fei::utils::cpu_time() - fei_creatematrix_start_time;
256 if (localProc == 0) {
257 std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
258 std::cout << "Create-Matrix cpu time: " << fei_creatematrix_time << std::endl;
259 }
260
261 double start_load_time = fei::utils::cpu_time();
262
263 fei::SharedPtr<fei::Vector> solnVec = factory->createVector(matrixGraph, true);
264 fei::SharedPtr<fei::Vector> rhsVec = factory->createVector(matrixGraph);
266 factory->createLinearSystem(matrixGraph);
267
268 CHK_ERR( linSys->parameters(params));
269
270 linSys->setMatrix(mat);
271 linSys->setSolutionVector(solnVec);
272 linSys->setRHS(rhsVec);
273
275 load_elem_data(matrixGraph.get(), mat.get(), rhsVec.get(), hexcube) );
276
277 //temporarily disable boundary-conditions
278 // CHK_ERR( load_BC_data(linSys, hexcube) );
279
280 if (!slave_constraints) {
282 load_constraints(linSys.get(), hexcube, firstLocalCRID) );
283 }
284
285 CHK_ERR( linSys->loadComplete() );
286
287 double fei_load_time = fei::utils::cpu_time() - start_load_time;
288
289 if (localProc == 0) {
290 //IOS macros are defined in fei_stdinc.h
291 std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
292 std::cout << "Coef. loading cpu time: " << fei_load_time << std::endl;
293 std::cout << "Total assembly wall time: "
294 << fei_init_time + fei_creatematrix_time + fei_load_time << std::endl;
295 }
296
297 //
298 //now the load phase is complete, so we're ready to launch the underlying
299 //solver and solve Ax=b
300 //
301
302 //First, check whether the params (set from the input .i file) specify
303 //a "Trilinos_Solver" string (which may be Amesos...). If not, it won't
304 //matter but if so, the factory may switch based on this value, when
305 //creating the solver object instance.
306
307 std::string solver_name_value;
308 params.getStringParamValue("Trilinos_Solver", solver_name_value);
309
310 const char* charptr_solvername =
311 solver_name_value.empty() ? NULL : solver_name_value.c_str();
312
313 fei::SharedPtr<fei::Solver> solver = factory->createSolver(charptr_solvername);
314
315 int status;
316 int itersTaken = 0;
317
318 if (localProc==0) std::cout << "solve..." << std::endl;
319 double start_solve_time = fei::utils::cpu_time();
320
321 int err = solver->solve(linSys.get(),
322 NULL, //preconditioningMatrix
323 params,
324 itersTaken,
325 status);
326
327 double solve_time = fei::utils::cpu_time()-start_solve_time;
328
329 if (err!=0 && localProc==0) {
330 std::cout << "solve returned err: " << err <<", status: "
331 << status << std::endl;
332 }
333
334 if (localProc==0) {
335 std::cout << " cpu-time in solve: " << solve_time << std::endl;
336 }
337
338 CHK_ERR( solnVec->scatterToOverlap() );
339
340 delete hexcubeptr;
341
342 double elapsed_cpu_time = fei::utils::cpu_time() - start_time;
343
344 //The following IOS_... macros are defined in base/fei_iostream.hpp
345 std::cout.setf(IOS_FIXED, IOS_FLOATFIELD);
346 if (localProc==0) {
347 std::cout << "Proc0 cpu times (seconds):" << std::endl
348 << " FEI initialize: " << fei_init_time << std::endl
349 << " FEI create-matrix: " << fei_creatematrix_time << std::endl
350 << " FEI load: " << fei_load_time << std::endl
351 << " solve: " << solve_time << std::endl
352 << "Total program time: " << elapsed_cpu_time << std::endl;
353
354 }
355
356
357#ifndef FEI_SER
358 MPI_Finalize();
359#endif
360
361 return(0);
362}
363
int init_elem_connectivities(FEI *fei, PoissonData &poissonData)
int load_elem_data(FEI *fei, PoissonData &poissonData)
int main(int argc, char **argv)
Definition: beam.cpp:85
@ OneD
Definition: HexBeam.hpp:26
virtual int numDofPerNode()
Definition: HexBeam.hpp:36
int numLocalDOF_
Definition: HexBeam.hpp:91
int totalNumElems_
Definition: HexBeam.hpp:75
virtual int getNumCRs()
Definition: HexBeam.hpp:63
int localNumElems_
Definition: HexBeam.hpp:77
int numGlobalDOF_
Definition: HexBeam.hpp:92
int getIntParamValue(const char *name, int &paramValue) const
int getStringParamValue(const char *name, std::string &paramValue) const
T * get() const
#define CHK_ERR(a)
#define IOS_FIXED
#define IOS_FLOATFIELD
#define MPI_COMM_WORLD
Definition: fei_mpi.h:58
#define MPI_Comm
Definition: fei_mpi.h:56
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
Definition: fei_utils.cpp:191
const char * version()
Definition: fei_utils.hpp:53
double cpu_time()
Definition: fei_utils.cpp:46
int initialize_mpi(int argc, char **argv, int &localProc, int &numProcs)
int get_filename_and_read_input(int argc, char **argv, MPI_Comm comm, int localProc, std::vector< std::string > &stdstrings)
fei::SharedPtr< fei::Factory > create_fei_Factory(MPI_Comm comm, const char *libraryName)
std::ostream & console_out()