MFEM  v4.6.0
Finite element discretization library
elasticity_operator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include "elasticity_operator.hpp"
14 
15 namespace mfem
16 {
17 
19  : Operator(), mesh_(mesh), order_(order), dim_(mesh_.SpaceDimension()),
20  vdim_(mesh_.SpaceDimension()), ne_(mesh_.GetNE()), h1_fec_(order_, dim_),
21  h1_fes_(&mesh_, &h1_fec_, vdim_, Ordering::byNODES)
22 {
23  this->height = h1_fes_.GetTrueVSize();
24  this->width = this->height;
25 
26  int global_tdof_size = h1_fes_.GlobalTrueVSize();
27  if (mesh.GetMyRank() == 0)
28  {
29  out << "#dofs: " << global_tdof_size << std::endl;
30  }
31 
35 
36  ir_ = const_cast<IntegrationRule *>(
38 
42  d1d_ = maps_->ndof;
43  q1d_ = maps_->nqpt;
44 
45  dX_ess_.UseDevice(true);
47 
48  X_el_.UseDevice(true);
50 
51  Y_el_.UseDevice(true);
53 
54  cstate_el_.UseDevice(true);
56 
57  X_local_.UseDevice(true);
59 
60  Y_local_.UseDevice(true);
62 
65 
66  if (use_cache_)
67  {
69  }
70 
72 }
73 
74 void ElasticityOperator::Mult(const Vector &X, Vector &Y) const
75 {
77 
78  // T-vector to L-vector
80  // L-vector to E-vector
82 
83  // Reset output vector
84  Y_el_ = 0.0;
85 
86  // Apply operator
89  X_el_, Y_el_);
90 
91  // E-vector to L-vector
93  // L-vector to T-vector
95 
96  // Set the residual at Dirichlet dofs on the T-vector to zero
98 }
99 
101 {
102  // invalidate cache
103  recompute_cache_ = true;
104 
107  return *gradient_;
108 }
109 
111 {
113 
114  // Column elimination for essential dofs
115  dX_ess_ = dX;
117 
118  // T-vector to L-vector
120  // L-vector to E-vector
122 
123  // Reset output vector
124  Y_el_ = 0.0;
125 
126  // Apply operator
130 
131  // E-vector to L-vector
133  // L-vector to T-vector
135 
136  // Re-assign the essential degrees of freedom on the final output vector.
137  {
138  const auto d_dX = dX.Read();
139  auto d_Y = Y.ReadWrite();
140  const auto d_ess_tdof_list = ess_tdof_list_.Read();
141  mfem::forall(ess_tdof_list_.Size(), [=] MFEM_HOST_DEVICE (int i)
142  {
143  d_Y[d_ess_tdof_list[i]] = d_dX[d_ess_tdof_list[i]];
144  });
145  }
146 
147  recompute_cache_ = false;
148 }
149 
151  Vector &K_diag_local,
152  Vector &K_diag) const
153 {
154  Ke_diag.SetSize(d1d_ * d1d_ * d1d_ * dim_ * ne_ * dim_);
155  K_diag_local.SetSize(h1_element_restriction_->Width() * dim_);
156  K_diag.SetSize(h1_prolongation_->Width() * dim_);
157 
160  geometric_factors_->detJ, cstate_el_, Ke_diag);
161 
162  // For each dimension, the H1 element restriction and H1 prolongation
163  // transpose actions are applied separately.
164  for (int i = 0; i < dim_; i++)
165  {
166  // Scalar component E-size
167  int sce_sz = d1d_ * d1d_ * d1d_ * dim_ * ne_;
168  // Scalar component L-size
169  int scl_sz = h1_element_restriction_->Width();
170 
171  Vector vin_local, vout_local;
172  vin_local.MakeRef(Ke_diag, i * sce_sz, sce_sz);
173  vout_local.MakeRef(K_diag_local, i * scl_sz, scl_sz);
174  h1_element_restriction_->MultTranspose(vin_local, vout_local);
175  vout_local.GetMemory().SyncAlias(K_diag_local.GetMemory(),
176  vout_local.Size());
177 
178  // Scalar component T-size
179  int sct_sz = h1_prolongation_->Width();
180  Vector vout;
181  vout.MakeRef(K_diag, i * sct_sz, sct_sz);
182  h1_prolongation_->MultTranspose(vout_local, vout);
183  vout.GetMemory().SyncAlias(K_diag.GetMemory(), vout.Size());
184  }
185 
186  // Each essential dof row and column are set to zero with it's diagonal entry
187  // set to 1, i.e. (Ke)_ii = 1.0.
189  int num_submats = h1_fes_.GetTrueVSize() / h1_fes_.GetVDim();
190  auto K_diag_submats = Reshape(K_diag.HostWrite(), num_submats, dim_, dim_);
191  for (int i = 0; i < ess_tdof_list_.Size(); i++)
192  {
193  int ess_idx = ess_tdof_list_[i];
194  int submat = ess_idx % num_submats;
195  int row = ess_idx / num_submats;
196  for (int j = 0; j < dim_; j++)
197  {
198  if (row == j)
199  {
200  K_diag_submats(submat, row, j) = 1.0;
201  }
202  else
203  {
204  K_diag_submats(submat, row, j) = 0.0;
205  K_diag_submats(submat, j, row) = 0.0;
206  }
207  }
208  }
209 }
210 
212 
213 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
std::function< void(const int, const Array< double > &, const Array< double > &, const Array< double > &, const Vector &, const Vector &, const Vector &, Vector &)> element_kernel_assemble_diagonal_wrapper
Wrapper for the assembly of the gradient on each diagonal element.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1152
const GeometricFactors * geometric_factors_
Vector X_el_
Input state E-vector.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:161
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:847
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:173
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:465
virtual void Mult(const Vector &U, Vector &Y) const override
Compute the residual Y = R(U) representing the elasticity equation with a material model chosen by ca...
ElasticityOperator(ParMesh &mesh, const int order)
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:169
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local, Vector &K_diag) const
Assemble the linearization of the residual K = dR(U)/dU.
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:228
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1302
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:365
std::function< void(const int, const Array< double > &, const Array< double > &, const Array< double > &, const Vector &, const Vector &, const Vector &, Vector &, const Vector &)> element_apply_gradient_kernel_wrapper
Wrapper for the application of the gradient of the residual.
Vector J
Jacobians of the element transformations at all quadrature points.
Definition: mesh.hpp:2224
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:86
Vector Y_el_
Output state E-Vector.
Vector detJ
Determinants of the Jacobians at all quadrature points.
Definition: mesh.hpp:2230
void SyncAlias(const Memory &base, int alias_size) const
Update the alias Memory *this to match the memory location (all valid locations) of its base Memory...
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
int q1d_
Number of quadrature points in 1D.
const Operator * h1_element_restriction_
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:311
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Vector X_local_
Input state L-vector.
Vector Y_local_
Output state L-vector.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
int GetMyRank() const
Definition: pmesh.hpp:353
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
const int ne_
Number of elements in the mesh (rank local)
void forall(int N, lambda &&body)
Definition: forall.hpp:742
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
int GetOrder(int i) const
Returns the polynomial degree of the i&#39;th finite element.
Definition: fespace.hpp:692
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:29
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:184
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
ElasticityGradientOperator is a wrapper class to pass ElasticityOperator::AssembleGradientDiagonal an...
int d1d_
Number of degrees of freedom in 1D.
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
std::function< void(const int, const Array< double > &, const Array< double > &, const Array< double > &, const Vector &, const Vector &, const Vector &, Vector &)> element_apply_kernel_wrapper
Wrapper for the application of the residual R(U).
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition: fe_base.hpp:205
Lexicographic ordering for tensor-product FiniteElements.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
Vector data type.
Definition: vector.hpp:58
void GradientMult(const Vector &dX, Vector &Y) const
Multiply the linearization of the residual R(U) wrt to the current state U by a perturbation dX...
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:581
ParFiniteElementSpace h1_fes_
H1 finite element space.
Abstract operator.
Definition: operator.hpp:24
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
Class for parallel meshes.
Definition: pmesh.hpp:32
Operator & GetGradient(const Vector &U) const override
Get the Gradient object.
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28