MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
elasticity_operator.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
14
15namespace 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
56
57 X_local_.UseDevice(true);
59
60 Y_local_.UseDevice(true);
62
65
66 if (use_cache_)
67 {
69 }
70
72}
73
74void 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
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 * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition array.hpp:321
int Size() const
Return the logical size of the array.
Definition array.hpp:144
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:317
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:210
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:161
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:189
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:174
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
ElasticityGradientOperator is a wrapper class to pass ElasticityOperator::AssembleGradientDiagonal an...
const Operator * h1_element_restriction_
int d1d_
Number of degrees of freedom in 1D.
Vector Y_local_
Output state L-vector.
Vector X_el_
Input state E-vector.
ElasticityOperator(ParMesh &mesh, const int order)
Vector Y_el_
Output state E-Vector.
std::function< void(const int, const Array< real_t > &, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, const Vector &, Vector &)> element_kernel_assemble_diagonal_wrapper
Wrapper for the assembly of the gradient on each diagonal element.
Vector X_local_
Input state L-vector.
const int ne_
Number of elements in the mesh (rank local)
int q1d_
Number of quadrature points in 1D.
std::function< void(const int, const Array< real_t > &, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, const Vector &, Vector &, const Vector &)> element_apply_gradient_kernel_wrapper
Wrapper for the application of the gradient of the residual.
std::function< void(const int, const Array< real_t > &, const Array< real_t > &, const Array< real_t > &, const Vector &, const Vector &, const Vector &, Vector &)> element_apply_kernel_wrapper
Wrapper for the application of the residual R(U).
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.
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...
Operator & GetGradient(const Vector &U) const override
Get the Gradient object.
const GeometricFactors * geometric_factors_
ParFiniteElementSpace h1_fes_
H1 finite element space.
void AssembleGradientDiagonal(Vector &Ke_diag, Vector &K_diag_local, Vector &K_diag) const
Assemble the linearization of the residual K = dR(U)/dU.
int GetOrder(int i) const
Returns the polynomial degree of the i'th finite element.
Definition fespace.hpp:696
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1336
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
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
Vector detJ
Determinants of the Jacobians at all quadrature points.
Definition mesh.hpp:2835
Vector J
Jacobians of the element transformations at all quadrature points.
Definition mesh.hpp:2829
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition intrules.cpp:86
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
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,...
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:856
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
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
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition fespace.hpp:30
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
const Operator * GetProlongationMatrix() const override
The returned Operator is owned by the FiniteElementSpace.
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
const FiniteElement * GetFE(int i) const override
Definition pfespace.cpp:534
Class for parallel meshes.
Definition pmesh.hpp:34
int GetMyRank() const
Definition pmesh.hpp:404
Vector data type.
Definition vector.hpp:80
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:474
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:490
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
Memory< real_t > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition vector.hpp:249
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:136
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:486
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
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
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
@ byNODES
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
@ LEXICOGRAPHIC
Lexicographic ordering for tensor-product FiniteElements.
void forall(int N, lambda &&body)
Definition forall.hpp:754
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486