Microsphere.hxx 13.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
/*!
//
// \file
//
//
// Created by Jérôme Diaz <jerome.diaz@inria.fr> on the Tue, 7 May 2019 16:35:18 +0100
// Copyright (c) Inria. All rights reserved.
//
//  Modified by Patrick Le Tallec to handle microsphere models
// Add microsphere contribution to free energy derivatives dW and d2W if  the local input fiber length are non zero.
// Only works in 3D.
// Data
// two family of fibers I4 and I6. For each family, need at node level
// Must also provide the integration rule on the sphere (number of points and angles) in radian respectively in angle phi and angle theta.
// \ingroup OperatorInstancesGroup
// \addtogroup OperatorInstancesGroup
// \{
*/


#ifndef MOREFEM_x_OPERATOR_INSTANCES_x_VARIATIONAL_OPERATOR_x_NONLINEAR_FORM_x_LOCAL_x_SECOND_PIOLA_KIRCHHOFF_STRESS_TENSOR_x_INTERNAL_VARIABLE_POLICY_x_MICROSPHERE_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_MICROSPHERE_HXX_


namespace MoReFEM
{


    namespace Advanced
    {


    namespace LocalVariationalOperatorNS
    {


        namespace SecondPiolaKirchhoffStressTensorNS
        {


            namespace InternalVariablePolicyNS
            {


                template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
                Microsphere<FiberIndexI4T, FiberIndexI6T>::Microsphere(const unsigned int mesh_dimension,
47 48 49 50
                                                                       const unsigned int Nnode_unknown,
                                                                       const unsigned int Nquad_point,
                                                                       const TimeManager& time_manager,
                                                                       input_internal_variable_policy_type* input_internal_variable_policy)
51 52 53 54 55 56 57 58 59 60 61
                : matrix_parent(),
                vector_parent(),
                time_manager_(time_manager),
                fibers_I4_(::MoReFEM::Internal::FiberNS::FiberListManager<ParameterNS::Type::vector>::GetInstance(__FILE__, __LINE__).GetFiberList(FiberIndexI4T)),
                fibers_I6_(::MoReFEM::Internal::FiberNS::FiberListManager<ParameterNS::Type::vector>::GetInstance(__FILE__, __LINE__).GetFiberList(FiberIndexI6T)),
                input_internal_variable_policy_(*input_internal_variable_policy)
                {
                    static_cast<void>(Nquad_point);
                    static_cast<void>(Nnode_unknown);
                    
                    unsigned int tauxtau_init(0u);
62
                    unsigned int tau_n_init(0u);
63 64
                    unsigned int tau_global_init(0u);

65 66
                    switch (mesh_dimension)
                    // Only the 3D case is relevant for this operator.
67 68 69
                    {
                        case 3u:
                            tauxtau_init = 6u;
70
                            tau_n_init = 3u;
71 72 73 74
                            tau_global_init = 3u;
                            break;
                        default:
                            assert(false);
75
							exit(EXIT_FAILURE);
76
                            break;
77
                    } // switch
78 79 80 81 82 83 84 85 86

                    this->matrix_parent::InitLocalMatrixStorage
                    ({{
                        {tauxtau_init, tauxtau_init}
                    }});

                    this->vector_parent::InitLocalVectorStorage
                    ({{
                        tauxtau_init,
87
                        tau_n_init,
88 89 90
                        tau_global_init
                    }});

91
                    // Hardcoded quadrature rule on the unit sphere.
92 93 94
                    const double pi = static_cast<double>(M_PI);
                    const auto factor_conv = pi / 180.;

95 96 97 98
					constexpr auto npquad = Utilities::ArraySize<decltype(lbdvp_)>::GetValue();

                    for (auto i = 0u; i < npquad; ++i)
                        lbdvp_[i] = pi * static_cast<double>(i + 1.) / static_cast<double>(npquad);
99

100
					constexpr auto ntquad = Utilities::ArraySize<decltype(lbdvt_)>::GetValue();
101

102 103
                    for (auto i = 0u; i < ntquad; ++i)
                        lbdvt_[i] = factor_conv * (70. + 10. * static_cast<double>(i));
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
                }


                template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
                template<unsigned int DimensionT>
                void Microsphere<FiberIndexI4T, FiberIndexI6T>
                ::ComputeWDerivates(const Advanced::LocalVariationalOperatorNS::InfosAtQuadPointNS::ForUnknownList& quad_pt_unknown_list_data,
                                    const GeometricElt& geom_elt,
                                    const Advanced::RefFEltInLocalOperator& ref_felt,
                                    const LocalVector& cauchy_green_tensor_value,
                                    const LocalMatrix& transposed_De,
                                    LocalVector& dW,
                                    LocalMatrix& d2W)
                {
                    static_cast<void>(ref_felt);
                    static_cast<void>(transposed_De);

                    const auto& quad_pt = quad_pt_unknown_list_data.GetQuadraturePoint();

123
                    // References to initiliazed vectors and matrices.
124 125 126 127 128 129
                    auto& tauxtau =
						this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::tauxtau)>();
                    auto& tau_n =
						this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::tau_n)>();
                    auto& tau_global =
						this->vector_parent::template GetLocalVector<EnumUnderlyingType(LocalVectorIndex::tau_global)>();
130

131
                    // Non linear stiffness contribution (second order derivative of the energetic potential).
132 133
                    auto& tauxtauxtauxtau =
						this->matrix_parent::template GetLocalMatrix<EnumUnderlyingType(LocalMatrixIndex::tauxtauxtauxtau)>();
134

135 136
                    // Get local values of fiber vector, fiber density, in plane fiber dispersion and out of plane fiber
                    // dispersion for families I4 and I6 from the InputPolicy. The fiber densities need to be prenormalized.
137 138
					decltype(auto) input_microsphere = GetInputMicrosphere();

139
                    const auto& fiber_I4 = fibers_I4_.GetValue(quad_pt, geom_elt);
140 141 142 143 144 145
                    const auto& density_I4 =
						input_microsphere.GetFiberStiffnessDensityI4().GetValue(quad_pt, geom_elt);
                    const auto& inplane_dispersion_I4 =
						input_microsphere.GetInPlaneFiberDispersionI4().GetValue(quad_pt, geom_elt);
                    const auto& transverse_dispersion_I4 =
						input_microsphere.GetOutOfPlaneFiberDispersionI4().GetValue(quad_pt, geom_elt);
146

147
                    const auto& fiber_I6 = fibers_I6_.GetValue(quad_pt, geom_elt);
148 149 150 151 152 153
                    const auto& density_I6 =
						input_microsphere.GetFiberStiffnessDensityI6().GetValue(quad_pt, geom_elt);
                    const auto& inplane_dispersion_I6 =
						input_microsphere.GetInPlaneFiberDispersionI6().GetValue(quad_pt, geom_elt);
                    const auto& transverse_dispersion_I6 =
						input_microsphere.GetOutOfPlaneFiberDispersionI6().GetValue(quad_pt, geom_elt);
154 155 156 157

                    // normalisation of fiber directions and calculation of the normal to the fiber plane
                    // fiber vectors are copied
					constexpr int norm2 = 2;
158 159 160 161
                    const auto norm_I4 = xt::linalg::norm(fiber_I4, norm2);
                    const auto norm_I6 = xt::linalg::norm(fiber_I6, norm2);
					auto tau_I4 = fiber_I4;
					auto tau_I6 = fiber_I6;
162

163 164 165 166 167 168 169 170 171 172 173 174
                    if (!(NumericNS::IsZero(norm_I4)) and !(NumericNS::IsZero(norm_I6)))
					{
						// local basis vectors tau_I4 & tau_I6
						tau_I4 /= norm_I4;
						tau_I6 /= norm_I6;
						
						// Vectorial product to construct third local basis vector tau_n (normal out of the plane).
						tau_n(0) = tau_I4(1) * tau_I6(2) - tau_I4(2) * tau_I6(1);
						tau_n(1) = tau_I4(2) * tau_I6(0) - tau_I4(0) * tau_I6(2);
						tau_n(2) = tau_I4(0) * tau_I6(1) - tau_I4(1) * tau_I6(0);
						
						const auto size = dW.shape(0);
175 176 177 178

						decltype(auto) quad_pt_list_along_plane = GetQuadraturePointsAlongPlane();
						decltype(auto) quad_pt_list_outside_plane = GetQuadraturePointsOutsidePlane();

179 180
						
						// Loop on spherical quadrature points.
181
						for (const auto value_quad_pt_p : quad_pt_list_along_plane)
182
						{
183
							for (const auto value_quad_pt_t : quad_pt_list_outside_plane)
184 185 186 187
							{
								double norm = 0.;
								
								// Local coordinates of the integration direction in the fiber frame of reference.
188 189 190
								const double xloc1 = std::cos(value_quad_pt_p) * std::sin(value_quad_pt_t);
								const double xloc2 = std::sin(value_quad_pt_p) * std::sin(value_quad_pt_t);
								const double xloc3 = std::cos(value_quad_pt_t);
191 192 193 194 195 196 197 198 199 200 201 202
								
								// Calculation of the local integration directions and diagonal values sigma_xx, sigma_yy and sigma_zz.
								for (auto i = 0u; i < DimensionT ; ++i)
								{
									tau_global(i) = tau_I4(i) * xloc1 + tau_I6(i) * xloc2 + tau_n(i) * xloc3;
									tauxtau(i) = NumericNS::Square(tau_global(i));
									norm += NumericNS::Square(tau_global(i));
								}
								
								// Extra diagonal values sigma_xy, sigma_yz and sigma_xz.
								for (auto i = DimensionT ; i < size ; ++i)
									tauxtau(i) = tau_global(i % DimensionT) * tau_global((i + 1) % DimensionT);
203

204 205
								// normalisation of initial fiberlength
								tauxtau /= norm;
206

GILLES Sebastien's avatar
GILLES Sebastien committed
207
								// calculation of the squ are elongation, of the fiber elongation, of the integration weights, of weigthed stiffness
208
                                double elongation_squared = xt::linalg::vdot(cauchy_green_tensor_value, tauxtau);
209

210 211
                                for (auto i = DimensionT; i < size ; ++i)
                                    elongation_squared += cauchy_green_tensor_value(i) * tauxtau(i);
212
								
213 214 215 216 217 218
								const double elongation = std::sqrt(elongation_squared);
								const double weight_I4 = std::exp(inplane_dispersion_I4 * std::cos(2.0 * value_quad_pt_p)
															- transverse_dispersion_I4 * std::cos(2.0 * value_quad_pt_t));
								const double weight_I6 = std::exp(- inplane_dispersion_I6 * std::cos(2.0 * value_quad_pt_p)
															- transverse_dispersion_I6 * std::cos(2.0 * value_quad_pt_t));
								const  double local_stiffness = weight_I4 * density_I4 + weight_I6 * density_I6;
219 220 221 222
								
								
								// Calculation of the fiber traction and of its derivative with respect to the squared elongation.
								// Piola = 2 * dW
223
								const double dW_contribution = local_stiffness * (elongation - 1.) / elongation;
224
								// Tangent = 4 * d2W
225
								const double d2W_contribution = local_stiffness / NumericNS::Cube(elongation);
226 227 228 229 230 231 232 233 234
								
								// update of gradient and of second derivative with respect to Cauchy Greene
								dW += dW_contribution * tauxtau;
								xt::noalias(tauxtauxtauxtau) = xt::linalg::outer(tauxtau, tauxtau);
								d2W += d2W_contribution * tauxtauxtauxtau;
							}
						}
					}
				}
235
				
236 237

                template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
238 239
                inline const ::MoReFEM::FiberList<ParameterNS::Type::vector>&
				Microsphere<FiberIndexI4T, FiberIndexI6T>::GetFibersI4() const noexcept
240 241 242 243 244
                {
                    return fibers_I4_;
                }

                template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
245 246
                inline const ::MoReFEM::FiberList<ParameterNS::Type::vector>&
				Microsphere<FiberIndexI4T, FiberIndexI6T>::GetFibersI6() const noexcept
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
                {
                    return fibers_I6_;
                }

                template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
                inline const TimeManager& Microsphere<FiberIndexI4T, FiberIndexI6T>::GetTimeManager() const noexcept
                {
                    return time_manager_;
                }


                template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
                inline const typename Microsphere<FiberIndexI4T, FiberIndexI6T>::input_internal_variable_policy_type&
                Microsphere<FiberIndexI4T, FiberIndexI6T>::GetInputMicrosphere() const noexcept
                {
                    return input_internal_variable_policy_;
                }


266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
				template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
                inline const std::array<double, 20ul>& Microsphere<FiberIndexI4T, FiberIndexI6T>
				::GetQuadraturePointsAlongPlane() const noexcept
				{
					return lbdvp_;
				}


				template <unsigned int FiberIndexI4T, unsigned int FiberIndexI6T>
				inline const std::array<double, 5ul>& Microsphere<FiberIndexI4T, FiberIndexI6T>
				::GetQuadraturePointsOutsidePlane() const noexcept
				{
					return lbdvt_;
				}


282
            } // namespace InternalVariablePolicyNS
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300


        } // namespace SecondPiolaKirchhoffStressTensorNS


    } // namespace LocalVariationalOperatorNS


    } // namespace Advanced


} // namespace MoReFEM


/// @} // addtogroup OperatorInstancesGroup


#endif // MOREFEM_x_OPERATOR_INSTANCES_x_VARIATIONAL_OPERATOR_x_NONLINEAR_FORM_x_LOCAL_x_SECOND_PIOLA_KIRCHHOFF_STRESS_TENSOR_x_INTERNAL_VARIABLE_POLICY_x_MICROSPHERE_HXX_