MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
generate_random_field.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details
11//
12// -------------------------------------------------------------------
13// Gaussian Random Fields of Matern Covariance for Imperfect Materials
14// -------------------------------------------------------------------
15//
16// See README.md for detailed description.
17//
18// Compile with: make generate_random_field
19//
20// Sample runs:
21// (Basic usage)
22// mpirun -np 4 generate_random_field
23//
24// (Generate 5 particles with random imperfections)
25// mpirun -np 4 generate_random_field -o 1 -r 3 -rp 3 -nu 2 -l1 0.015 -l2 0.015 -l3 0.015 -s 0.01 -t 0.08 -n 5 -pl2 3 -top 0 -rs
26//
27// (Generate an Octet-Truss with random imperfections)
28// mpirun -np 4 generate_random_field -o 1 -r 3 -rp 3 -nu 2 -l1 0.02 -l2 0.02 -l3 0.02 -s 0.01 -t 0.08 -top 1 -rs
29//
30// (Generate an Octet-Truss with random imperfections following a uniform distribution)
31// mpirun -np 4 generate_random_field -o 1 -r 3 -rp 3 -nu 2 -l1 0.02 -l2 0.02 -l3 0.02 -umin 0.01 -umax 0.05 -t 0.08 -top 1 -urf -rs
32//
33// (2D random field with anisotropy)
34// mpirun -np 4 generate_random_field -o 1 -r 3 -rp 3 -nu 4 -l1 0.09 -l2 0.03 -l3 0.05 -s 0.01 -t 0.08 -top 1 -no-rs -m ../../data/ref-square.mesh
35
36#include <math.h>
37#include <fstream>
38#include <iostream>
39#include <string>
40#include "mfem.hpp"
41
42#include "material_metrics.hpp"
43#include "spde_solver.hpp"
44#include "transformation.hpp"
45#include "util.hpp"
46#include "visualizer.hpp"
47
48using namespace std;
49using namespace mfem;
50
52
53int main(int argc, char *argv[])
54{
55#ifdef MFEM_USE_SINGLE
56 cout << "This miniapp is not supported in single precision.\n\n";
57 return MFEM_SKIP_RETURN_VALUE;
58#endif
59
60 // 0. Initialize MPI.
61 Mpi::Init(argc, argv);
63
64 // 1. Parse command-line options.
65 const char *mesh_file = "../../data/ref-cube.mesh";
66 int order = 1;
67 int num_refs = 3;
68 int num_parallel_refs = 3;
69 int number_of_particles = 3;
70 int topological_support = TopologicalSupport::kOctetTruss;
71 real_t nu = 2.0;
72 real_t tau = 0.08;
73 real_t l1 = 0.02;
74 real_t l2 = 0.02;
75 real_t l3 = 0.02;
76 real_t e1 = 0;
77 real_t e2 = 0;
78 real_t e3 = 0;
79 real_t pl1 = 1.0;
80 real_t pl2 = 1.0;
81 real_t pl3 = 1.0;
82 real_t uniform_min = 0.0;
83 real_t uniform_max = 1.0;
84 real_t offset = 0.0;
85 real_t scale = 0.01;
86 real_t level_set_threshold = 0.0;
87 bool paraview_export = true;
88 bool glvis_export = true;
89 bool uniform_rf = false;
90 bool random_seed = true;
91 bool compute_boundary_integrals = false;
92
93 OptionsParser args(argc, argv);
94 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
95 args.AddOption(&order, "-o", "--order",
96 "Finite element order (polynomial degree) or -1 for"
97 " isoparametric space.");
98 args.AddOption(&num_refs, "-r", "--refs", "Number of uniform refinements");
99 args.AddOption(&num_parallel_refs, "-rp", "--refs-parallel",
100 "Number of uniform refinements");
101 args.AddOption(&topological_support, "-top", "--topology",
102 "Topological support. 0 particles, 1 octet-truss");
103 args.AddOption(&nu, "-nu", "--nu", "Fractional exponent nu (smoothness)");
104 args.AddOption(&tau, "-t", "--tau", "Parameter for topology generation");
105 args.AddOption(&l1, "-l1", "--l1",
106 "First component of diagonal core of theta");
107 args.AddOption(&l2, "-l2", "--l2",
108 "Second component of diagonal core of theta");
109 args.AddOption(&l3, "-l3", "--l3",
110 "Third component of diagonal core of theta");
111 args.AddOption(&e1, "-e1", "--e1", "First euler angle for rotation of theta");
112 args.AddOption(&e2, "-e2", "--e2",
113 "Second euler angle for rotation of theta");
114 args.AddOption(&e3, "-e3", "--e3", "Third euler angle for rotation of theta");
115 args.AddOption(&pl1, "-pl1", "--pl1", "Length scale 1 of particles");
116 args.AddOption(&pl2, "-pl2", "--pl2", "Length scale 2 of particles");
117 args.AddOption(&pl3, "-pl3", "--pl3", "Length scale 3 of particles");
118 args.AddOption(&uniform_min, "-umin", "--uniform-min",
119 "Minimum value of uniform distribution");
120 args.AddOption(&uniform_max, "-umax", "--uniform-max",
121 "Maximum value of uniform distribution");
122 args.AddOption(&offset, "-off", "--offset",
123 "Offset for random field u(x) -> u(x) + a");
124 args.AddOption(&scale, "-s", "--scale",
125 "Scale for random field u(x) -> a * u(x)");
126 args.AddOption(&level_set_threshold, "-lst", "--level-set-threshold",
127 "Level set threshold");
128 args.AddOption(&number_of_particles, "-n", "--number-of-particles",
129 "Number of particles");
130 args.AddOption(&paraview_export, "-pvis", "--paraview-visualization",
131 "-no-pvis", "--no-paraview-visualization",
132 "Enable or disable ParaView visualization.");
133 args.AddOption(&glvis_export, "-vis", "--visualization", "-no-vis",
134 "--no-visualization",
135 "Enable or disable GLVis visualization.");
136 args.AddOption(&uniform_rf, "-urf", "--uniform-rf", "-no-urf",
137 "--no-uniform-rf",
138 "Enable or disable the transformation of GRF to URF.");
139 args.AddOption(&random_seed, "-rs", "--random-seed", "-no-rs",
140 "--no-random-seed", "Enable or disable random seed.");
141 args.AddOption(&compute_boundary_integrals, "-cbi",
142 "--compute-boundary-integrals", "-no-cbi",
143 "--no-compute-boundary-integrals",
144 "Enable or disable computation of boundary integrals.");
145 args.Parse();
146 if (!args.Good())
147 {
148 args.PrintUsage(cout);
149 return 1;
150 }
151 if (Mpi::Root())
152 {
153 args.PrintOptions(cout);
154 }
155
156 // 2. Read the mesh from the given mesh file.
157 Mesh mesh(mesh_file, 1, 1);
158 int dim = mesh.Dimension();
159 bool is_3d = (dim == 3);
160
161 // 3. Refine the mesh to increase the resolution.
162 for (int i = 0; i < num_refs; i++)
163 {
164 mesh.UniformRefinement();
165 }
166 ParMesh pmesh(MPI_COMM_WORLD, mesh);
167 mesh.Clear();
168 for (int i = 0; i < num_parallel_refs; i++)
169 {
170 pmesh.UniformRefinement();
171 }
172
173 // 4. Define a finite element space on the mesh.
174 H1_FECollection fec(order, dim);
175 ParFiniteElementSpace fespace(&pmesh, &fec);
176 HYPRE_BigInt size = fespace.GlobalTrueVSize();
177 if (Mpi::Root())
178 {
179 const Array<int> boundary(pmesh.bdr_attributes);
180 cout << "Number of finite element unknowns: " << size << "\n";
181 cout << "Boundary attributes: ";
182 boundary.Print(cout, 6);
183 }
184
185 // ========================================================================
186 // II. Generate topological support
187 // ========================================================================
188
189 ParGridFunction v(&fespace);
190 v = 0.0;
191 MaterialTopology *mdm = nullptr;
192
193 // II.1 Define the metric for the topological support.
194 if (is_3d)
195 {
196 if (topological_support == TopologicalSupport::kOctetTruss)
197 {
198 mdm = new OctetTrussTopology();
199 }
200 else if (topological_support == TopologicalSupport::kParticles)
201 {
202 // Create the same random particles on all processors.
203 std::vector<real_t> random_positions(3 * number_of_particles);
204 std::vector<real_t> random_rotations(9 * number_of_particles);
205 if (Mpi::Root())
206 {
207 // Generate random positions and rotations. We generate them on the root
208 // process and then broadcast them to all processes because we need the
209 // same random positions and rotations on all processes.
210 FillWithRandomNumbers(random_positions, 0.2, 0.8);
211 FillWithRandomRotations(random_rotations);
212 }
213
214 // Broadcast the random positions and rotations to all processes.
215 MPI_Bcast(random_positions.data(), 3 * number_of_particles,
216 MPITypeMap<real_t>::mpi_type, 0, MPI_COMM_WORLD);
217 MPI_Bcast(random_rotations.data(), 9 * number_of_particles,
218 MPITypeMap<real_t>::mpi_type, 0, MPI_COMM_WORLD);
219
220 mdm = new ParticleTopology(pl1, pl2, pl3, random_positions,
221 random_rotations);
222 }
223 else
224 {
225 if (Mpi::Root())
226 {
227 mfem::out << "Error: Selected topological support not valid."
228 << std::endl;
229 }
230 return 1;
231 }
232
233 // II.2 Define lambda to wrap the call to the distance metric.
234 auto topo = [&mdm, &tau](const Vector &x)
235 {
236 return (tau - mdm->ComputeMetric(x));
237 };
238
239 // II.3 Create a GridFunction for the topological support.
240 FunctionCoefficient topo_coeff(topo);
241 v.ProjectCoefficient(topo_coeff);
242 }
243
244 // ========================================================================
245 // III. Generate random imperfections via fractional PDE
246 // ========================================================================
247
248 /// III.1 Define the fractional PDE solution
249 ParGridFunction u(&fespace);
250 u = 0.0;
251
252 // III.2 Define the boundary conditions.
254 if (Mpi::Root())
255 {
256 bc.PrintInfo();
257 bc.VerifyDefinedBoundaries(pmesh);
258 }
259
260 // III.3 Solve the SPDE problem
261 spde::SPDESolver solver(nu, bc, &fespace, l1, l2, l3, e1, e2,
262 e3);
263 const int seed = (random_seed) ? 0 : std::numeric_limits<int>::max();
264 solver.SetupRandomFieldGenerator(seed);
265 solver.GenerateRandomField(u);
266
267 /// III.4 Verify boundary conditions
268 if (compute_boundary_integrals)
269 {
271 }
272
273 // ========================================================================
274 // III. Combine topological support and random field
275 // ========================================================================
276
277 if (uniform_rf)
278 {
279 /// Transform the random field to a uniform random field.
280 spde::UniformGRFTransformer transformation(uniform_min, uniform_max);
281 transformation.Transform(u);
282 }
283 if (scale != 1.0)
284 {
285 /// Scale the random field.
287 transformation.Transform(u);
288 }
289 if (offset != 0.0)
290 {
291 /// Add an offset to the random field.
293 transformation.Transform(u);
294 }
295 ParGridFunction w(&fespace); // Noisy material field.
296 w = 0.0;
297 w += u;
298 w += v;
299 ParGridFunction level_set(w); // Level set field.
300 {
301 spde::LevelSetTransformer transformation(level_set_threshold);
302 transformation.Transform(level_set);
303 }
304
305 // ========================================================================
306 // IV. Export visualization to ParaView and GLVis
307 // ========================================================================
308
309 spde::Visualizer vis(pmesh, order, u, v, w, level_set, is_3d);
310 if (paraview_export)
311 {
312 vis.ExportToParaView();
313 }
314 if (glvis_export)
315 {
316 vis.SendToGLVis();
317 }
318
319 delete mdm;
320 return 0;
321}
void Print(std::ostream &out=mfem::out, int width=4) const
Prints array to stream with width elements per row.
Definition array.cpp:23
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Virtual class to define the interface for defining the material topology.
virtual real_t ComputeMetric(const Vector &x)=0
Compute the metric rho describing the material topology.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel meshes.
Definition pmesh.hpp:34
Class that implements the particle topology.
Vector data type.
Definition vector.hpp:80
Level Set Transformer, 1 for u(x) >= threshold, 0 otherwise.
Adds an constant offset to a grid function, i.e. u(x) = u(x) + offset.
void GenerateRandomField(ParGridFunction &x)
void SetupRandomFieldGenerator(int seed=0)
Transforms a grid function by scaling it by a constant factor.
void SendToGLVis() const
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
void transformation(const Vector &p, Vector &v)
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void FillWithRandomRotations(std::vector< real_t > &x)
Definition util.cpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
void FillWithRandomNumbers(std::vector< real_t > &x, real_t a, real_t b)
Fills the vector x with random numbers between a and b.
Definition util.cpp:18
float real_t
Definition config.hpp:43
Helper struct to convert a C++ type to an MPI type.
void PrintInfo(std::ostream &os=mfem::out) const
Print the information specifying the boundary conditions.
void ComputeBoundaryError(const ParGridFunction &solution)
void VerifyDefinedBoundaries(const Mesh &mesh) const