AnalyticalPrestress.hxx 11.1 KB
Newer Older
GILLES Sebastien's avatar
GILLES Sebastien committed
1 2 3 4 5 6 7 8 9 10 11 12 13
/*!
//
// \file
//
//
// Created by Gautier Bureau <gautier.bureau@inria.fr> on the Thu, 14 Jan 2016 12:00:52 +0100
// Copyright (c) Inria. All rights reserved.
//
// \ingroup OperatorInstancesGroup
// \addtogroup OperatorInstancesGroup
// \{
*/

14

15 16
#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_VARIATIONAL_OPERATOR_x_NONLINEAR_FORM_x_LOCAL_x_SECOND_PIOLA_KIRCHHOFF_STRESS_TENSOR_x_INTERNAL_VARIABLE_POLICY_x_ANALYTICAL_PRESTRESS_HXX_
# define MOREFEM_x_OPERATOR_INSTANCES_x_VARIATIONAL_OPERATOR_x_NONLINEAR_FORM_x_LOCAL_x_SECOND_PIOLA_KIRCHHOFF_STRESS_TENSOR_x_INTERNAL_VARIABLE_POLICY_x_ANALYTICAL_PRESTRESS_HXX_
17 18


19
namespace MoReFEM
20
{
21 22


23 24
    namespace Advanced
    {
25 26


27
    namespace LocalVariationalOperatorNS
28
    {
29 30


31
        namespace SecondPiolaKirchhoffStressTensorNS
32
        {
33 34


35
            namespace InternalVariablePolicyNS
36
            {
37 38


39
                template <unsigned int FiberIndexT>
40
                AnalyticalPrestress<FiberIndexT>::AnalyticalPrestress(const unsigned int mesh_dimension,
41 42 43
                                                                      const unsigned int Nnode_unknown,
                                                                      const unsigned int Nquad_point,
                                                                      const TimeManager& time_manager,
44
                                                                      input_internal_variable_policy_type* input_internal_variable_policy)
45
                : vector_parent(),
46
                fibers_(::MoReFEM::Internal::FiberNS::FiberListManager<ParameterNS::Type::vector>::GetInstance(__FILE__, __LINE__).GetFiberList(FiberIndexT)),
47
                time_manager_(time_manager),
48
                input_internal_variable_policy_(*input_internal_variable_policy)
49
                {
50
                    static_cast<void>(Nquad_point);
51

52
                    unsigned int tauxtau_init(0u);
53

54 55 56 57 58 59 60 61 62 63 64 65
                    switch(mesh_dimension)
                    {
                        case 2u:
                            tauxtau_init = 3u;
                            break;
                        case 3u:
                            tauxtau_init = 6u;
                            break;
                        default:
                            assert(false);
                            break;
                    }
66

67 68
                    this->vector_parent::InitLocalVectorStorage
                    ({{
69
                        tauxtau_init
70
                    }});
71

72 73
                    local_electrical_activation_previous_time_.resize(Nnode_unknown);
                    local_electrical_activation_at_time_.resize(Nnode_unknown);
74
                }
75 76


77 78
                template <unsigned int FiberIndexT>
                template<unsigned int DimensionT>
79
                void AnalyticalPrestress<FiberIndexT>
80
                ::ComputeWDerivates(const Advanced::LocalVariationalOperatorNS::InfosAtQuadPointNS::ForUnknownList& quad_pt_unknown_list_data,
81 82
                                    const GeometricElt& geom_elt,
                                    const Advanced::RefFEltInLocalOperator& ref_felt,
83
                                    const LocalVector& cauchy_green_tensor_value,
84 85 86
                                    const LocalMatrix& transposed_De,
                                    LocalVector& dW,
                                    LocalMatrix& d2W)
87 88
                {
                    static_cast<void>(d2W);
89
                    static_cast<void>(ref_felt);
90
                    static_cast<void>(cauchy_green_tensor_value);
91
                    static_cast<void>(transposed_De);
92

93
                    const auto run_case = GetTimeManager().GetStaticOrDynamic();
94

95
                    if (run_case == StaticOrDynamic::dynamic_)
96
                    {
97 98
                        const auto& quad_pt = quad_pt_unknown_list_data.GetQuadraturePoint();
                        const auto& phi_ref = quad_pt_unknown_list_data.GetRefFEltPhi();
99

100 101 102
                        auto& fibers = GetFibers();
                        auto& tauxtau = this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::tauxtau)>();
                        const double dt = GetTimeManager().GetTimeStep();
103

104
                        const auto& tau_interpolate = fibers.GetValue(quad_pt, geom_elt);
105

GILLES Sebastien's avatar
GILLES Sebastien committed
106
                        const auto dimension = static_cast<std::size_t>(DimensionT);
107

108
                        double norm = 0.;
109

GILLES Sebastien's avatar
GILLES Sebastien committed
110
                        for (auto component = 0ul; component < dimension ; ++component)
111
                            norm += NumericNS::Square(tau_interpolate(component));
112

GILLES Sebastien's avatar
GILLES Sebastien committed
113
                        const auto size = dW.shape(0);
114

115
                        // Diagonal values sigma_xx, sigma_yy and sigma_zz.
GILLES Sebastien's avatar
GILLES Sebastien committed
116
                        for (auto i = 0ul ; i < dimension ; ++i)
Gautier Bureau's avatar
Gautier Bureau committed
117
                            tauxtau(i) = NumericNS::Square(tau_interpolate(i));
118

119
                        // Extra diagonal values sigma_xy, sigma_yz and sigma_xz.
GILLES Sebastien's avatar
GILLES Sebastien committed
120
                        for (auto i = dimension ; i < size ; ++i)
121 122
                            tauxtau(i) = tau_interpolate(i % dimension) * tau_interpolate((i + 1) % dimension);

123 124
                        if (!(NumericNS::IsZero(norm)))
                        {
125
                            tauxtau /= norm;
126

Gautier Bureau's avatar
Gautier Bureau committed
127
                            double new_sigma_c = 0.;
128

Gautier Bureau's avatar
Gautier Bureau committed
129
                            if (do_update_sigma_c_)
130
                            {
Gautier Bureau's avatar
Gautier Bureau committed
131 132
                                const auto& u0 = GetLocalElectricalActivationPreviousTime();
                                const auto& u1 = GetLocalElectricalActivationAtTime();
133

Gautier Bureau's avatar
Gautier Bureau committed
134 135
                                double u0_at_quad = 0.;
                                double u1_at_quad = 0.;
136

137
                                const double sigma_0 = GetInputActiveStress().GetContractility().GetValue(quad_pt, geom_elt);
138

Gautier Bureau's avatar
Gautier Bureau committed
139
                                const int Nnode = static_cast<int>(geom_elt.Nvertex());
140

Gautier Bureau's avatar
Gautier Bureau committed
141
                                // This loop interpolates u at the current quad point.
142
                                for (auto node_index = 0; node_index < Nnode; ++node_index)
Gautier Bureau's avatar
Gautier Bureau committed
143 144 145 146
                                {
                                    u0_at_quad += phi_ref(node_index) * u0[static_cast<unsigned int>(node_index)];
                                    u1_at_quad += phi_ref(node_index) * u1[static_cast<unsigned int>(node_index)];
                                }
147

Gautier Bureau's avatar
Gautier Bureau committed
148
                                const double u0_u1_at_quad = u0_at_quad + u1_at_quad;
149

Gautier Bureau's avatar
Gautier Bureau committed
150 151 152 153 154 155
                                new_sigma_c = GetNonCstSigmaC().UpdateAndGetValue(quad_pt, geom_elt,
                                                                                  [u0_u1_at_quad, sigma_0, dt](double& sigma_c)
                                                                                  {
                                                                                      sigma_c = (0.5 * sigma_0 * NumericNS::AbsPlus(u0_u1_at_quad) + (1. / dt - 0.5 * std::fabs(u0_u1_at_quad)) * sigma_c) / (1. / dt + 0.5 * std::fabs(u0_u1_at_quad));
                                                                                  }
                                                                                  );
156
                            }
Gautier Bureau's avatar
Gautier Bureau committed
157 158
                            else
                                new_sigma_c = GetNonCstSigmaC().GetValue(quad_pt, geom_elt);
159

160
                            xt::noalias(dW) += new_sigma_c * tauxtau;
161
                        }
162 163
                    }
                }
164 165


166
                template <unsigned int FiberIndexT>
167
                inline const ::MoReFEM::FiberList<ParameterNS::Type::vector>& AnalyticalPrestress<FiberIndexT>::GetFibers() const noexcept
168 169 170
                {
                    return fibers_;
                }
171 172


173
                template <unsigned int FiberIndexT>
174
                inline const TimeManager& AnalyticalPrestress<FiberIndexT>::GetTimeManager() const noexcept
175 176 177
                {
                    return time_manager_;
                }
178 179


180
                template <unsigned int FiberIndexT>
181 182
                inline typename AnalyticalPrestress<FiberIndexT>::ScalarParameterAtQuadPt&
                AnalyticalPrestress<FiberIndexT>::GetNonCstSigmaC() noexcept
183
                {
184
                    return const_cast<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>&>(GetSigmaC());
185
                }
186 187


188
                template <unsigned int FiberIndexT>
189 190
                inline const typename AnalyticalPrestress<FiberIndexT>::ScalarParameterAtQuadPt&
                AnalyticalPrestress<FiberIndexT>::GetSigmaC() const noexcept
191 192 193 194
                {
                    assert(!(!sigma_c_));
                    return *sigma_c_;
                }
195 196


197
                template <unsigned int FiberIndexT>
198 199
                inline void AnalyticalPrestress<FiberIndexT>
                ::SetSigmaC(AnalyticalPrestress<FiberIndexT>::ScalarParameterAtQuadPt* sigma_c)
200 201 202
                {
                    sigma_c_ = sigma_c;
                }
203 204


205 206
                template <unsigned int FiberIndexT>
                inline const std::vector<double>&
207
                AnalyticalPrestress<FiberIndexT>::GetLocalElectricalActivationPreviousTime() const noexcept
208
                {
209
                    return local_electrical_activation_previous_time_;
210
                }
211 212


213 214
                template <unsigned int FiberIndexT>
                inline std::vector<double>&
215
                AnalyticalPrestress<FiberIndexT>::GetNonCstLocalElectricalActivationPreviousTime() noexcept
216
                {
217
                    return const_cast<std::vector<double>&>(GetLocalElectricalActivationPreviousTime());
218
                }
219 220


221 222
                template <unsigned int FiberIndexT>
                inline const std::vector<double>&
223
                AnalyticalPrestress<FiberIndexT>::GetLocalElectricalActivationAtTime() const noexcept
224
                {
225
                    return local_electrical_activation_at_time_;
226
                }
227 228


229 230
                template <unsigned int FiberIndexT>
                inline std::vector<double>&
231
                AnalyticalPrestress<FiberIndexT>::GetNonCstLocalElectricalActivationAtTime() noexcept
232
                {
233
                    return const_cast<std::vector<double>&>(GetLocalElectricalActivationAtTime());
234
                }
235 236


237
                template <unsigned int FiberIndexT>
238
                void AnalyticalPrestress<FiberIndexT>::SetDoUpdateSigmaC(const bool update)
239
                {
240
                    do_update_sigma_c_ = update;
241
                }
242 243


244
                template <unsigned int FiberIndexT>
245
                inline const typename AnalyticalPrestress<FiberIndexT>::input_internal_variable_policy_type&
246
                AnalyticalPrestress<FiberIndexT>::GetInputActiveStress() const noexcept
247
                {
248
                    return input_internal_variable_policy_;
249
                }
250 251


252
            } // namespace InternalVariablePolicyNS
253 254


255
        } // namespace SecondPiolaKirchhoffStressTensorNS
256 257


258
    } // namespace LocalVariationalOperatorNS
259 260


261
    } // namespace Advanced
262 263


264
} // namespace MoReFEM
265 266


267 268 269
/// @} // addtogroup OperatorInstancesGroup


270
#endif // MOREFEM_x_OPERATOR_INSTANCES_x_VARIATIONAL_OPERATOR_x_NONLINEAR_FORM_x_LOCAL_x_SECOND_PIOLA_KIRCHHOFF_STRESS_TENSOR_x_INTERNAL_VARIABLE_POLICY_x_ANALYTICAL_PRESTRESS_HXX_