forked from libMesh/libmesh
-
Notifications
You must be signed in to change notification settings - Fork 0
/
transient_ex2.C
662 lines (533 loc) · 23.5 KB
/
transient_ex2.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
/* The libMesh Finite Element Library. */
/* Copyright (C) 2003 Benjamin S. Kirk */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */
/* License as published by the Free Software Foundation; either */
/* version 2.1 of the License, or (at your option) any later version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */
/* Lesser General Public License for more details. */
/* You should have received a copy of the GNU Lesser General Public */
/* License along with this library; if not, write to the Free Software */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
// <h1>Transient Example 2 - The Newmark System and the Wave Equation</h1>
//
// This is the eighth example program. It builds on
// the previous example programs. It introduces the
// NewmarkSystem class. In this example the wave equation
// is solved using the time integration scheme provided
// by the NewmarkSystem class.
//
// This example comes with a cylindrical mesh given in the
// universal file pipe-mesh.unv.
// The mesh contains HEX8 and PRISM6 elements.
// C++ include files that we need
#include <iostream>
#include <fstream>
#include <algorithm>
#include <stdio.h>
#include <math.h>
// Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/gmv_io.h"
#include "libmesh/vtk_io.h"
#include "libmesh/newmark_system.h"
#include "libmesh/equation_systems.h"
// Define the Finite Element object.
#include "libmesh/fe.h"
// Define Gauss quadrature rules.
#include "libmesh/quadrature_gauss.h"
// Define useful datatypes for finite element
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
// Define matrix and vector data types for the global
// equation system. These are base classes,
// from which specific implementations, like
// the PETSc or LASPACK implementations, are derived.
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
// Define the DofMap, which handles degree of freedom
// indexing.
#include "libmesh/dof_map.h"
// The definition of a vertex associated with a Mesh.
#include "libmesh/node.h"
// The definition of a geometric element
#include "libmesh/elem.h"
// Defines the MeshData class, which allows you to store
// data about the mesh when reading in files, etc.
#include "libmesh/mesh_data.h"
// Bring in everything from the libMesh namespace
using namespace libMesh;
// Function prototype. This is the function that will assemble
// the linear system for our problem, governed by the linear
// wave equation.
void assemble_wave(EquationSystems& es,
const std::string& system_name);
// Function Prototype. This function will be used to apply the
// initial conditions.
void apply_initial(EquationSystems& es,
const std::string& system_name);
// Function Prototype. This function imposes
// Dirichlet Boundary conditions via the penalty
// method after the system is assembled.
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name);
// The main program
int main (int argc, char** argv)
{
// Initialize libraries, like in example 2.
LibMeshInit init (argc, argv);
// Check for proper usage.
if (argc < 2)
{
if (libMesh::processor_id() == 0)
std::cerr << "Usage: " << argv[0] << " [meshfile]"
<< std::endl;
libmesh_error();
}
// Tell the user what we are doing.
else
{
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
}
// LasPack solvers don't work so well for this example
// (not sure why), and Trilinos matrices don't work at all.
// Print a warning to the user if PETSc is not in use.
if (libMesh::default_solver_package() == LASPACK_SOLVERS)
{
std::cout << "WARNING! It appears you are using the\n"
<< "LasPack solvers. This example may not converge\n"
<< "using LasPack, but should work OK with PETSc.\n"
<< "http://www.mcs.anl.gov/petsc/\n"
<< std::endl;
}
else if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
{
std::cout << "WARNING! It appears you are using the\n"
<< "Trilinos solvers. The current libMesh Epetra\n"
<< "interface does not allow sparse matrix addition,\n"
<< "as is needed in this problem. We recommend\n"
<< "using PETSc: http://www.mcs.anl.gov/petsc/\n"
<< std::endl;
return 0;
}
// Get the name of the mesh file
// from the command line.
std::string mesh_file = argv[1];
std::cout << "Mesh file is: " << mesh_file << std::endl;
// Skip this 3D example if libMesh was compiled as 1D or 2D-only.
libmesh_example_assert(3 <= LIBMESH_DIM, "3D support");
// Create a mesh.
// This example directly references all mesh nodes and is
// incompatible with ParallelMesh use.
//
// Create a SerialMesh object, with dimension to be overridden
// later, distributed across the default MPI communicator.
SerialMesh mesh(init.communicator());
MeshData mesh_data(mesh);
// Read the meshfile specified in the command line or
// use the internal mesh generator to create a uniform
// grid on an elongated cube.
mesh.read(mesh_file, &mesh_data);
// mesh.build_cube (10, 10, 40,
// -1., 1.,
// -1., 1.,
// 0., 4.,
// HEX8);
// Print information about the mesh to the screen.
mesh.print_info();
// The node that should be monitored.
const unsigned int result_node = 274;
// Time stepping issues
//
// Note that the total current time is stored as a parameter
// in the \pEquationSystems object.
//
// the time step size
const Real delta_t = .0000625;
// The number of time steps.
unsigned int n_time_steps = 300;
// Create an equation systems object.
EquationSystems equation_systems (mesh);
// Declare the system and its variables.
// Create a NewmarkSystem named "Wave"
equation_systems.add_system<NewmarkSystem> ("Wave");
// Use a handy reference to this system
NewmarkSystem & t_system = equation_systems.get_system<NewmarkSystem> ("Wave");
// Add the variable "p" to "Wave". "p"
// will be approximated using first-order approximation.
t_system.add_variable("p", FIRST);
// Give the system a pointer to the matrix assembly
// function and the initial condition function defined
// below.
t_system.attach_assemble_function (assemble_wave);
t_system.attach_init_function (apply_initial);
// Set the time step size, and optionally the
// Newmark parameters, so that \p NewmarkSystem can
// compute integration constants. Here we simply use
// pass only the time step and use default values
// for \p alpha=.25 and \p delta=.5.
t_system.set_newmark_parameters(delta_t);
// Set the speed of sound and fluid density
// as \p EquationSystems parameter,
// so that \p assemble_wave() can access it.
equation_systems.parameters.set<Real>("speed") = 1000.;
equation_systems.parameters.set<Real>("fluid density") = 1000.;
// Start time integration from t=0
t_system.time = 0.;
// Initialize the data structures for the equation system.
equation_systems.init();
// Prints information about the system to the screen.
equation_systems.print_info();
// A file to store the results at certain nodes.
std::ofstream res_out("pressure_node.res");
// get the dof_numbers for the nodes that
// should be monitored.
const unsigned int res_node_no = result_node;
const Node& res_node = mesh.node(res_node_no-1);
unsigned int dof_no = res_node.dof_number(0,0,0);
// Assemble the time independent system matrices and rhs.
// This function will also compute the effective system matrix
// K~=K+a_0*M+a_1*C and apply user specified initial
// conditions.
t_system.assemble();
// Now solve for each time step.
// For convenience, use a local buffer of the
// current time. But once this time is updated,
// also update the \p EquationSystems parameter
// Start with t_time = 0 and write a short header
// to the nodal result file
res_out << "# pressure at node " << res_node_no << "\n"
<< "# time\tpressure\n"
<< t_system.time << "\t" << 0 << std::endl;
for (unsigned int time_step=0; time_step<n_time_steps; time_step++)
{
// Update the time. Both here and in the
// \p EquationSystems object
t_system.time += delta_t;
// Update the rhs.
t_system.update_rhs();
// Impose essential boundary conditions.
// Not that since the matrix is only assembled once,
// the penalty parameter should be added to the matrix
// only in the first time step. The applied
// boundary conditions may be time-dependent and hence
// the rhs vector is considered in each time step.
if (time_step == 0)
{
// The local function \p fill_dirichlet_bc()
// may also set Dirichlet boundary conditions for the
// matrix. When you set the flag as shown below,
// the flag will return true. If you want it to return
// false, simply do not set it.
equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = true;
fill_dirichlet_bc(equation_systems, "Wave");
// unset the flag, so that it returns false
equation_systems.parameters.set<bool>("Newmark set BC for Matrix") = false;
}
else
fill_dirichlet_bc(equation_systems, "Wave");
// Solve the system "Wave".
t_system.solve();
// After solving the system, write the solution
// to a GMV-formatted plot file.
// Do only for a few time steps.
if (time_step == 30 || time_step == 60 ||
time_step == 90 || time_step == 120 )
{
char buf[14];
if (!libMesh::on_command_line("--vtk")){
sprintf (buf, "out.%03d.gmv", time_step);
GMVIO(mesh).write_equation_systems (buf,equation_systems);
}else{
#ifdef LIBMESH_HAVE_VTK
// VTK viewers are generally not happy with two dots in a filename
sprintf (buf, "out_%03d.exd", time_step);
VTKIO(mesh).write_equation_systems (buf,equation_systems);
#endif // #ifdef LIBMESH_HAVE_VTK
}
}
// Update the p, v and a.
t_system.update_u_v_a();
// dof_no may not be local in parallel runs, so we may need a
// global displacement vector
NumericVector<Number> &displacement
= t_system.get_vector("displacement");
std::vector<Number> global_displacement(displacement.size());
displacement.localize(global_displacement);
// Write nodal results to file. The results can then
// be viewed with e.g. gnuplot (run gnuplot and type
// 'plot "pressure_node.res" with lines' in the command line)
res_out << t_system.time << "\t"
<< global_displacement[dof_no]
<< std::endl;
}
// All done.
return 0;
}
// This function assembles the system matrix and right-hand-side
// for our wave equation.
void assemble_wave(EquationSystems& es,
const std::string& system_name)
{
// It is a good idea to make sure we are assembling
// the proper system.
libmesh_assert_equal_to (system_name, "Wave");
// Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
// The dimension that we are running.
const unsigned int dim = mesh.mesh_dimension();
// Copy the speed of sound and fluid density
// to a local variable.
const Real speed = es.parameters.get<Real>("speed");
const Real rho = es.parameters.get<Real>("fluid density");
// Get a reference to our system, as before.
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
// Get a constant reference to the Finite Element type
// for the first (and only) variable in the system.
FEType fe_type = t_system.get_dof_map().variable_type(0);
// In here, we will add the element matrices to the
// @e additional matrices "stiffness_mass" and "damping"
// and the additional vector "force", not to the members
// "matrix" and "rhs". Therefore, get writable
// references to them.
SparseMatrix<Number>& stiffness = t_system.get_matrix("stiffness");
SparseMatrix<Number>& damping = t_system.get_matrix("damping");
SparseMatrix<Number>& mass = t_system.get_matrix("mass");
NumericVector<Number>& force = t_system.get_vector("force");
// Some solver packages (PETSc) are especially picky about
// allocating sparsity structure and truly assigning values
// to this structure. Namely, matrix additions, as performed
// later, exhibit acceptable performance only for identical
// sparsity structures. Therefore, explicitly zero the
// values in the collective matrix, so that matrix additions
// encounter identical sparsity structures.
SparseMatrix<Number>& matrix = *t_system.matrix;
DenseMatrix<Number> zero_matrix;
// Build a Finite Element object of the specified type. Since the
// \p FEBase::build() member dynamically creates memory we will
// store the object as an \p AutoPtr<FEBase>. This can be thought
// of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
// A 2nd order Gauss quadrature rule for numerical integration.
QGauss qrule (dim, SECOND);
// Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
// The element Jacobian * quadrature weight at each integration point.
const std::vector<Real>& JxW = fe->get_JxW();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
// The element shape function gradients evaluated at the quadrature
// points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
// A reference to the \p DofMap object for this system. The \p DofMap
// object handles the index translation from node and element numbers
// to degree of freedom numbers.
const DofMap& dof_map = t_system.get_dof_map();
// The element mass, damping and stiffness matrices
// and the element contribution to the rhs.
DenseMatrix<Number> Ke, Ce, Me;
DenseVector<Number> Fe;
// This vector will hold the degree of freedom indices for
// the element. These define where in the global system
// the element degrees of freedom get mapped.
std::vector<dof_id_type> dof_indices;
// Now we will loop over all the elements in the mesh.
// We will compute the element matrix and right-hand-side
// contribution.
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
// Store a pointer to the element we are currently
// working on. This allows for nicer syntax later.
const Elem* elem = *el;
// Get the degree of freedom indices for the
// current element. These define where in the global
// matrix and right-hand-side this element will
// contribute to.
dof_map.dof_indices (elem, dof_indices);
// Compute the element-specific data for the current
// element. This involves computing the location of the
// quadrature points (q_point) and the shape functions
// (phi, dphi) for the current element.
fe->reinit (elem);
// Zero the element matrices and rhs before
// summing them. We use the resize member here because
// the number of degrees of freedom might have changed from
// the last element. Note that this will be the case if the
// element type is different (i.e. the last element was HEX8
// and now have a PRISM6).
{
const unsigned int n_dof_indices = dof_indices.size();
Ke.resize (n_dof_indices, n_dof_indices);
Ce.resize (n_dof_indices, n_dof_indices);
Me.resize (n_dof_indices, n_dof_indices);
zero_matrix.resize (n_dof_indices, n_dof_indices);
Fe.resize (n_dof_indices);
}
// Now loop over the quadrature points. This handles
// the numeric integration.
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
// Now we will build the element matrix. This involves
// a double loop to integrate the test funcions (i) against
// the trial functions (j).
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp]
*1./(speed*speed);
} // end of the matrix summation loop
} // end of quadrature point loop
// Now compute the contribution to the element matrix and the
// right-hand-side vector if the current element lies on the
// boundary.
{
// In this example no natural boundary conditions will
// be considered. The code is left here so it can easily
// be extended.
//
// don't do this for any side
for (unsigned int side=0; side<elem->n_sides(); side++)
if (!true)
// if (elem->neighbor(side) == NULL)
{
// Declare a special finite element object for
// boundary integration.
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
// Boundary integration requires one quadraure rule,
// with dimensionality one less than the dimensionality
// of the element.
QGauss qface(dim-1, SECOND);
// Tell the finte element object to use our
// quadrature rule.
fe_face->attach_quadrature_rule (&qface);
// The value of the shape functions at the quadrature
// points.
const std::vector<std::vector<Real> >& phi_face = fe_face->get_phi();
// The Jacobian * Quadrature Weight at the quadrature
// points on the face.
const std::vector<Real>& JxW_face = fe_face->get_JxW();
// Compute the shape function values on the element
// face.
fe_face->reinit(elem, side);
// Here we consider a normal acceleration acc_n=1 applied to
// the whole boundary of our mesh.
const Real acc_n_value = 1.0;
// Loop over the face quadrature points for integration.
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
// Right-hand-side contribution due to prescribed
// normal acceleration.
for (unsigned int i=0; i<phi_face.size(); i++)
{
Fe(i) += acc_n_value*rho
*phi_face[i][qp]*JxW_face[qp];
}
} // end face quadrature point loop
} // end if (elem->neighbor(side) == NULL)
// In this example the Dirichlet boundary conditions will be
// imposed via panalty method after the
// system is assembled.
} // end boundary condition section
// If this assembly program were to be used on an adaptive mesh,
// we would have to apply any hanging node constraint equations
// by uncommenting the following lines:
// std::vector<unsigned int> dof_indicesC = dof_indices;
// std::vector<unsigned int> dof_indicesM = dof_indices;
// dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
// dof_map.constrain_element_matrix (Ce, dof_indicesC);
// dof_map.constrain_element_matrix (Me, dof_indicesM);
// Finally, simply add the contributions to the additional
// matrices and vector.
stiffness.add_matrix (Ke, dof_indices);
damping.add_matrix (Ce, dof_indices);
mass.add_matrix (Me, dof_indices);
force.add_vector (Fe, dof_indices);
// For the overall matrix, explicitly zero the entries where
// we added values in the other ones, so that we have
// identical sparsity footprints.
matrix.add_matrix(zero_matrix, dof_indices);
} // end of element loop
// All done!
return;
}
// This function applies the initial conditions
void apply_initial(EquationSystems& es,
const std::string& system_name)
{
// Get a reference to our system, as before
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
// Numeric vectors for the pressure, velocity and acceleration
// values.
NumericVector<Number>& pres_vec = t_system.get_vector("displacement");
NumericVector<Number>& vel_vec = t_system.get_vector("velocity");
NumericVector<Number>& acc_vec = t_system.get_vector("acceleration");
// Assume our fluid to be at rest, which would
// also be the default conditions in class NewmarkSystem,
// but let us do it explicetly here.
pres_vec.zero();
vel_vec.zero();
acc_vec.zero();
}
// This function applies the Dirichlet boundary conditions
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name)
{
// It is a good idea to make sure we are assembling
// the proper system.
libmesh_assert_equal_to (system_name, "Wave");
// Get a reference to our system, as before.
NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
// Get writable references to the overall matrix and vector.
SparseMatrix<Number>& matrix = *t_system.matrix;
NumericVector<Number>& rhs = *t_system.rhs;
// Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
// Get \p libMesh's \f$ \pi \f$
const Real pi = libMesh::pi;
// Ask the \p EquationSystems flag whether
// we should do this also for the matrix
const bool do_for_matrix =
es.parameters.get<bool>("Newmark set BC for Matrix");
// Number of nodes in the mesh.
unsigned int n_nodes = mesh.n_nodes();
for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
{
// Get a reference to the current node.
const Node& curr_node = mesh.node(n_cnt);
// Check if Dirichlet BCs should be applied to this node.
// Use the \p TOLERANCE from \p mesh_common.h as tolerance.
// Here a pressure value is applied if the z-coord.
// is equal to 4, which corresponds to one end of the
// pipe-mesh in this directory.
const Real z_coo = 4.;
if (fabs(curr_node(2)-z_coo) < TOLERANCE)
{
// The global number of the respective degree of freedom.
unsigned int dn = curr_node.dof_number(0,0,0);
// The penalty parameter.
const Real penalty = 1.e10;
// Here we apply sinusoidal pressure values for 0<t<0.002
// at one end of the pipe-mesh.
Real p_value;
if (t_system.time < .002 )
p_value = sin(2*pi*t_system.time/.002);
else
p_value = .0;
// Now add the contributions to the matrix and the rhs.
rhs.add(dn, p_value*penalty);
// Add the panalty parameter to the global matrix
// if desired.
if (do_for_matrix)
matrix.add(dn, dn, penalty);
}
} // loop n_cnt
}