67 : lambda(lambda_),
mu(mu_),
u(NULL), si(0), sj(0) { }
70 void SetComponent(
int i,
int j) { si = i; sj = j; }
76 class VisMan :
public iostream
84 int win_x, win_y, win_w, win_h;
85 int win_stride_x, win_stride_y, win_nx;
90 void CloseConnection();
91 void PositionWindow();
99 ostream &
operator<<(ostream &v,
void (*
f)(VisMan&));
101 int main(
int argc,
char *argv[])
103 Mpi::Init(argc, argv);
107 const char *mesh_file =
"../data/beam-tri.mesh";
108 int ser_ref_levels = -1;
109 int par_ref_levels = 1;
113 bool amg_elast =
false;
114 bool visualization = 1;
117 args.
AddOption(&mesh_file,
"-m",
"--mesh",
118 "Mesh file to use.");
119 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
120 "Number of times to refine the mesh uniformly before parallel"
121 " partitioning, -1 for auto.");
122 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
123 "Number of times to refine the mesh uniformly after parallel"
126 "Finite element order (polynomial degree).");
128 "One of the two DG penalty parameters, typically +1/-1."
129 " See the documentation of class DGElasticityIntegrator.");
131 "One of the two DG penalty parameters, should be positive."
132 " Negative values are replaced with (order+1)^2.");
133 args.
AddOption(&amg_elast,
"-elast",
"--amg-for-elasticity",
"-sys",
135 "Use the special AMG elasticity solver (GM/LN approaches), "
136 "or standard AMG for systems (unknown approach).");
137 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
138 "--no-visualization",
139 "Enable or disable GLVis visualization.");
148 kappa = (order+1)*(order+1);
153 Mesh mesh(mesh_file, 1, 1);
160 cerr <<
"\nInput mesh should have at least two materials and "
161 <<
"two boundary attributes! (See schematic in ex17p.cpp)\n"
168 if (ser_ref_levels < 0)
172 for (
int l = 0; l < ser_ref_levels; l++)
180 ParMesh pmesh(MPI_COMM_WORLD, mesh);
183 for (
int l = 0; l < par_ref_levels; l++)
197 cout <<
"Number of finite element unknowns: " << glob_size
198 <<
"\nAssembling: " << flush;
239 if (Mpi::Root()) { cout <<
"r.h.s. ... " << flush; }
242 init_x, lambda_c, mu_c, alpha, kappa), dir_bdr);
260 if (Mpi::Root()) { cout <<
"matrix ... " << flush; }
266 if (Mpi::Root()) { cout <<
"done." << endl; }
271 const double rtol = 1e-6;
301 if (visualization) { reference_nodes = *pmesh.
GetNodes(); }
310 ostringstream mesh_name, sol_name;
311 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << Mpi::WorldRank();
312 sol_name <<
"sol." << setfill(
'0') << setw(6) << Mpi::WorldRank();
314 ofstream mesh_ofs(mesh_name.str().c_str());
315 mesh_ofs.precision(8);
318 ofstream sol_ofs(sol_name.str().c_str());
319 sol_ofs.precision(8);
328 VisMan vis(vishost, visport);
329 const char *glvis_keys = (dim < 3) ?
"Rjlc" :
"c";
335 <<
"solution\n" << pmesh << x << flush
336 <<
"keys " << glvis_keys << endl
337 <<
"window_title 'Deformed configuration'" << endl
338 <<
"plot_caption 'Backward displacement'" << endl
342 const char *c =
"xyz";
345 StressCoefficient stress_c(lambda_c, mu_c);
346 *pmesh.
GetNodes() = reference_nodes;
348 stress_c.SetDisplacement(x);
349 for (
int si = 0; si <
dim; si++)
351 for (
int sj = si; sj <
dim; sj++)
353 stress_c.SetComponent(si, sj);
356 MPI_Barrier(MPI_COMM_WORLD);
360 <<
"solution\n" << pmesh << stress << flush
361 <<
"keys " << glvis_keys << endl
362 <<
"window_title |Stress " << c[si] << c[sj] <<
'|' << endl
375 u(u.
Size()-1) = -0.2*x(0);
382 MFEM_ASSERT(u != NULL,
"displacement field is not set");
384 double L = lambda.Eval(T, ip);
385 double M =
mu.Eval(T, ip);
386 u->GetVectorGradient(T, grad);
389 double div_u = grad.Trace();
390 return L*div_u + 2*M*grad(si,si);
394 return M*(grad(si,sj) + grad(sj,si));
401 host(vishost), port(visport), sid(0)
407 win_stride_x = win_w;
408 win_stride_y = win_h + 20;
412 void VisMan::NewWindow()
416 iostream::rdbuf(sock[sid]->rdbuf());
419 void VisMan::CloseConnection()
421 if (sid < sock.Size())
429 void VisMan::PositionWindow()
431 *
this <<
"window_geometry "
432 << win_x + win_stride_x*(sid%win_nx) <<
' '
433 << win_y + win_stride_y*(sid/win_nx) <<
' '
434 << win_w <<
' ' << win_h << endl;
439 for (
int i = sock.Size()-1; i >= 0; i--)
445 ostream &
operator<<(ostream &v,
void (*
f)(VisMan&))
447 VisMan *vp =
dynamic_cast<VisMan*
>(&v);
448 if (vp) { (*f)(*vp); }
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Evaluate the coefficient in the element described by T at the point ip at time t. ...
Conjugate gradient method.
MPI_Comm GetComm() const
MPI communicator.
Class for grid function - Vector with associated FE space.
HYPRE_BigInt GlobalTrueVSize() const
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
int GetNE() const
Returns number of elements.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Abstract parallel finite element space.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void position_window(VisMan &v)
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
The BoomerAMG solver in hypre.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void new_window(VisMan &v)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void SetMaxIter(int max_it)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
void SetElasticityOptions(ParFiniteElementSpace *fespace)
void SetNodalFESpace(FiniteElementSpace *nfes)
void PrintUsage(std::ostream &out) const
Print the usage message.
A general vector function coefficient.
void close_connection(VisMan &v)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
Abstract base class for iterative solver.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
NURBSExtension * NURBSext
Optional NURBS mesh extension.
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
Class for integration point with weight.
void InitDisplacement(const Vector &x, Vector &u)
void PrintOptions(std::ostream &out) const
Print the options.
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
void Clear()
Clear the contents of the Mesh.
void GetNodes(Vector &node_coord) const
void SetSystemsOptions(int dim, bool order_bynodes=false)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
double u(const Vector &xvec)
Class for parallel grid function.
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
void Neg()
(*this) = -(*this)
bool Good() const
Return true if the command line options were parsed successfully.