diff --git a/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.C b/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.C index 3f844fb21..5e95a8081 100644 --- a/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.C +++ b/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.C @@ -39,7 +39,8 @@ #include "classfactory.h" #include "mathfem.h" #include "tensor/tensor3.h" - +#include "domain.h" +#include "function.h" namespace oofem { @@ -56,21 +57,14 @@ HardMagneticMooneyRivlinCompressibleMaterial :: give_FirstPKStressVector_Magneti Tensor1_3d H(vH), Q, B; Tensor2_3d F(vF), P; //Q = (H+M) - Q(i_3) = H(i_3) + M(i_3); + double m_level = this->giveDomain()->giveFunction( m_ltf )->evaluateAtTime( tStep->giveIntrinsicTime() ); + Q(i_3) = H(i_3) + m_level * M(i_3); // compute the first Piola-Kirchhoff auto [J, cofF] = F.compute_determinant_and_cofactor(); // P( i_3, j_3 ) = C1 * this->compute_dI1_Cdev_dF(F)(i_3, j_3) + C2 * this->compute_dI2_Cdev_dF(F)(i_3, j_3) + this->compute_dVolumetricEnergy_dF(F)(i_3, j_3) - - ( mu_0 / J ) * F.compute_tensor_cross_product()( i_3, j_3, m_3, n_3 ) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 )); - //+ (mu_0 / (2*J*J)) * cofF(m_3, n_3) * Q(n_3) * cofF(m_3, k_3) * Q(k_3) * cofF(i_3, j_3); - - /* - Tensor2_3d TT, TL; - TT( i_3, j_3 ) = - ( mu_0 / J ) * F.compute_tensor_cross_product()( i_3, j_3, m_3, n_3 ) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 )); - Tensor3_3d eps; - eps.be_Levi_Civita(); - TL( i_3, j_3 ) = - ( mu_0 / J ) * eps(i_3, k_3, m_3) * F( k_3, l_3) * eps(j_3, l_3, n_3) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 )); - */ + - ( mu_0 / J ) * F.compute_tensor_cross_product()( i_3, j_3, m_3, n_3 ) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 )) + + (mu_0 / (2*J*J)) * cofF(m_3, n_3) * Q(n_3) * cofF(m_3, k_3) * Q(k_3) * cofF(i_3, j_3); B( i_3 ) = mu_0 / J * cofF( j_3, i_3 ) * cofF( j_3, k_3 ) * Q( k_3 ); auto vP = P.to_voigt_form(); @@ -110,24 +104,51 @@ HardMagneticMooneyRivlinCompressibleMaterial :: giveConstitutiveMatrices_dPdF_dB // //Q = (H+M) - Q(i_3) = H(i_3) + M(i_3); + double m_level = this->giveDomain()->giveFunction( m_ltf )->evaluateAtTime( tStep->giveIntrinsicTime() ); + Q(i_3) = H(i_3) + m_level * M(i_3); GQQ( i_3, j_3 ) = G( i_3, k_3 ) * Q( k_3 ) * Q( j_3 ); auto GQQcross = GQQ.compute_tensor_cross_product(); - // - + // Tensor4_3d dPdF; Tensor3_3d dPdH; Tensor3_3d dBdF; Tensor2_3d dBdH; dPdF(i_3, j_3, k_3, l_3) = C1 * this->compute_d2I1_Cdev_dF2(F)(i_3, j_3, k_3, l_3) - + C2 * this->compute_d2I2_Cdev_dF2(F)(i_3, j_3, k_3, l_3) - + this->compute_d2VolumetricEnergy_dF2(F)(i_3, j_3, k_3, l_3) + + C2 * this->compute_d2I2_Cdev_dF2(F)(i_3, j_3, k_3, l_3) + + this->compute_d2VolumetricEnergy_dF2(F)(i_3, j_3, k_3, l_3) - mu_0 / J * (GQQcross(i_3, j_3, k_3, l_3) + ( Fcross( i_3, j_3, m_3, n_3 ) * Q(n_3) ) * ( Fcross( m_3, q_3, k_3, l_3 ) * Q( q_3 ) )) - + mu_0/J/J * Fcross( i_3, j_3, m_3, n_3 ) * GQQ(m_3, n_3) * G(k_3, l_3) ; - // + ( mu_0 / ( 2 * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q( k_3 ) * Fcross( i_3, j_3, k_3, l_3 ) ) - // - ( mu_0 / ( J * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q(k_3) * G(i_3, j_3) * G(k_3, l_3) - // + ( mu_0 / ( J * J ) ) * G(i_3, j_3) * (Fcross(k_3, l_3, m_3, n_3) * G(m_3,q_3) * Q(q_3) * Q(n_3))); + + mu_0/J/J * Fcross( i_3, j_3, m_3, n_3 ) * GQQ(m_3, n_3) * G(k_3, l_3) + + ( mu_0 / ( 2 * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q( k_3 ) * Fcross( i_3, j_3, k_3, l_3 ) ) + - ( mu_0 / ( J * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q(k_3) * G(i_3, j_3) * G(k_3, l_3)) + + ( mu_0 / ( J * J ) ) * G(i_3, j_3) * (Fcross(k_3, l_3, m_3, n_3) * G(m_3,q_3) * Q(q_3) * Q(n_3)); // + + Tensor4_3d A1, A2, A2t, A3; + //A1(i_3, j_3, k_3, l_3) = + ( mu_0 / ( 2 * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q( k_3 ) * Fcross( i_3, j_3, k_3, l_3 ) ); + //A1(i_3, j_3, k_3, l_3) = - ( mu_0 / ( J * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q(k_3) * G(i_3, j_3) * G(k_3, l_3); + // A1(i_3, j_3, k_3, l_3) = + ( mu_0 / ( J * J ) ) * G(i_3, j_3) * (Fcross(k_3, l_3, m_3, n_3) * G(m_3,q_3) * Q(q_3) * Q(n_3))); + + /* + A1(i_3, j_3, k_3, l_3) = - mu_0 / J * GQQcross(i_3, j_3, k_3, l_3); + Tensor3_3d FcQ1, FcQ2, FcQ3, FcQ4; + FcQ1(i_3, j_3, m_3) = Fcross( i_3, j_3, m_3, n_3 ) * Q(n_3); + FcQ2(m_3, k_3, l_3) = Fcross( m_3, q_3, k_3, l_3 ) * Q( q_3 ); + FcQ3(m_3, q_3, l_3) = Q( k_3 ) * Fcross( m_3, q_3, k_3, l_3 ); + FcQ4(m_3, q_3, k_3) = Q( l_3 ) * Fcross( m_3, q_3, k_3, l_3 ); + auto Fc = FloatMatrix(Fcross.to_voigt_form()); + + auto FcQ1m = FloatMatrix(FcQ1.to_voigt_form_3x9()); + auto FcQ2m = FloatMatrix(FcQ2.to_voigt_form_3x9()); + auto FcQ3m = FloatMatrix(FcQ3.to_voigt_form_3x9()); + auto FcQ4m = FloatMatrix(FcQ4.to_voigt_form_3x9()); + A2t(i_3, j_3, k_3, l_3) = FcQ1(i_3, j_3, m_3) * FcQ2(m_3, k_3, l_3); + A2(i_3, j_3, k_3, l_3) =- mu_0 / J * ( Fcross( i_3, j_3, m_3, n_3 ) * Q(n_3) ) * ( Fcross( m_3, q_3, k_3, l_3 ) * Q( q_3 ) ); + A3(i_3, j_3, k_3, l_3) = + mu_0/J/J * Fcross( i_3, j_3, m_3, n_3 ) * GQQ(m_3, n_3) * G(k_3, l_3) ; + // + auto A1v = FloatMatrix(A1.to_voigt_form()); + auto A2v = FloatMatrix(A2.to_voigt_form()); + auto A3v = FloatMatrix(A3.to_voigt_form()); + */ dBdH( i_3, k_3 ) = -mu_0 / J * G( j_3, i_3 ) * G( j_3, k_3 ); // dBdF( i_3, k_3, l_3 ) = mu_0 / J * (G( j_3, m_3 ) * Q( m_3 )* Fcross( j_3, i_3, k_3, l_3 )) + G( j_3, i_3 ) * (Q( m_3 ) * Fcross( j_3, m_3, k_3, l_3 )) - mu_0 / ( J * J ) * ( G( j_3, i_3 ) * G( j_3, m_3 ) * Q( m_3 ) * G( k_3, l_3 ) ); @@ -139,7 +160,7 @@ HardMagneticMooneyRivlinCompressibleMaterial :: giveConstitutiveMatrices_dPdF_dB auto [K1, K2, K3, K4] = this->compute_stiff_num(vF, vH, gp, tStep); - + auto D1 = FloatMatrix(dPdF.to_voigt_form()); auto D2 = FloatMatrix(dPdH.to_voigt_form_9x3()); auto D3 = FloatMatrix(dBdF.to_voigt_form_3x9()); @@ -223,11 +244,40 @@ HardMagneticMooneyRivlinCompressibleMaterial::initializeFrom(InputRecord &ir) IR_GIVE_FIELD(ir, C2, _IFT_HardMagneticMooneyRivlinCompressibleMaterial_c2); FloatArray m; IR_GIVE_FIELD(ir, m, _IFT_HardMagneticMooneyRivlinCompressibleMaterial_magnetization); + IR_GIVE_FIELD(ir, m_ltf, _IFT_HardMagneticMooneyRivlinCompressibleMaterial_m_ltf); if(m.giveSize() != 3){ OOFEM_ERROR("Magnetization has to be vector of size 3"); } else { M = Tensor1_3d(FloatArrayF<3>(m)); } } + + + + +int +HardMagneticMooneyRivlinCompressibleMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) +{ + + MagnetoElasticMaterialStatus *status = static_cast< MagnetoElasticMaterialStatus * >( this->giveStatus(gp) ); + if ( type == IST_DeformationGradientTensor ) { + answer = status->giveFVector(); + return 1; + } else if ( type == IST_FirstPKStressTensor ) { + answer = status->givePVector(); + return 1; + } else if ( type == IST_MagneticInductionVector ) { + answer = status->giveBVector(); + return 1; + } else if ( type == IST_MagneticFieldVector ) { + answer = status->giveHVector(); + return 1; + } else { + return Material::giveIPValue(answer, gp, type, tStep); + } +} + + + } // end namespace oofem diff --git a/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.h b/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.h index ba42f6734..07165fedf 100644 --- a/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.h +++ b/src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.h @@ -44,6 +44,7 @@ //@{ #define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_Name "hardmagneticmooneyrivlincompressiblemat" #define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_magnetization "magnetization" +#define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_m_ltf "m_ltf" #define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_c1 "c1" #define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_c2 "c2" @@ -65,6 +66,7 @@ class HardMagneticMooneyRivlinCompressibleMaterial : public MagnetoElasticMateri double C1; double C2; double mu_0 = 1.25663706143e-6; + int m_ltf = 0; @@ -86,6 +88,10 @@ class HardMagneticMooneyRivlinCompressibleMaterial : public MagnetoElasticMateri const char *giveInputRecordName() const override { return _IFT_HardMagneticMooneyRivlinCompressibleMaterial_Name; } const char *giveClassName() const override { return "HardMagneticMooneyRivlinCompressibleMaterial"; } + + int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep); + + }; } // end namespace oofem #endif diff --git a/src/mpm/magnetoelasticelement.C b/src/mpm/magnetoelasticelement.C index 0c7015570..7ac579e75 100644 --- a/src/mpm/magnetoelasticelement.C +++ b/src/mpm/magnetoelasticelement.C @@ -42,6 +42,7 @@ #include "intarray.h" #include "classfactory.h" #include "gaussintegrationrule.h" +#include "crosssection.h" #include "fei2dquadquad.h" #include "fei2dquadlin.h" @@ -95,7 +96,8 @@ class MagnetoElasticElement : public MPElement { } else { OOFEM_ERROR("Unknown characteristic vector type"); } - } + } + }; @@ -185,5 +187,83 @@ REGISTER_Element(MagnetoElasticQuad_ql) +/** + * @brief 2D Magneto elatic element with linear interpolation of displacements and magnetic potential + * + */ +class MagnetoElasticQuad_ll : public MagnetoElasticElement { + protected: + const static FEInterpolation & phiInterpol; + const static FEInterpolation & uInterpol; + const static Variable& phi; + const static Variable& u; + + public: + MagnetoElasticQuad_ll(int n, Domain* d): MagnetoElasticElement(n,d) + { + numberOfDofMans = 4; + numberOfGaussPoints = 4; + this->computeGaussPoints(); + } + + int getNumberOfSurfaceDOFs() const override {return 0;} + int getNumberOfEdgeDOFs() const override {return 0;} + void getSurfaceLocalCodeNumbers(oofem::IntArray&, oofem::Variable::VariableQuantity) const override {;} + + void getDofManLocalCodeNumbers (IntArray& answer, const Variable::VariableQuantity q, int num ) const override { + /* dof ordering: u1 v1 phi1 u2 v2 phi2 u3 v3 phi3 u4 v4 phi4 u5 v5 u6 v6 */ + if (q == Variable::VariableQuantity::Displacement) { + //answer={1,2, 4,5, 7,8, 10,11, 13,14 15,16, 17,18 19,20}; + int o = (num-1)*3+1; + answer={o, o+1}; + } else if (q == Variable::VariableQuantity::MagneticPotential) { + answer={num*3}; + } + } + void getInternalDofManLocalCodeNumbers (IntArray& answer, const Variable::VariableQuantity q, int num ) const override { + answer={}; + } + + void giveDofManDofIDMask(int inode, IntArray &answer) const override { + answer = {1,2,55}; + } + int giveNumberOfDofs() override { return 12; } + const char *giveInputRecordName() const override {return "MagnetoElasticQuad_ll";} + const char *giveClassName() const override { return "MagnetoElasticQuad_ll"; } + + + const FEInterpolation& getGeometryInterpolation() const override {return this->uInterpol;} + + Element_Geometry_Type giveGeometryType() const override { + return EGT_quad_1; + } + void getEdgeLocalCodeNumbers(IntArray& answer, const Variable::VariableQuantity q) const override {} + + +private: + virtual int giveNumberOfUDofs() const override {return 8;} + virtual int giveNumberOfPhiDofs() const override {return 4;} + virtual const Variable& getU() const override {return u;} + virtual const Variable& getPhi() const override {return phi;} + void computeGaussPoints() override { + if ( integrationRulesArray.size() == 0 ) { + integrationRulesArray.resize( 1 ); + integrationRulesArray [ 0 ] = std::make_unique(1, this); + integrationRulesArray [ 0 ]->SetUpPointsOnSquare(numberOfGaussPoints, _PlaneStrain); + } + } +}; + + const FEInterpolation & MagnetoElasticQuad_ll :: uInterpol = FEI2dQuadLin(1,2); + const FEInterpolation & MagnetoElasticQuad_ll :: phiInterpol = FEI2dQuadLin(1,2); +const Variable& MagnetoElasticQuad_ll::phi = Variable( MagnetoElasticQuad_ll::phiInterpol, Variable::VariableQuantity::MagneticPotential, Variable::VariableType::scalar, 1, NULL, {55}); +const Variable& MagnetoElasticQuad_ll::u = Variable(MagnetoElasticQuad_ll::uInterpol, Variable::VariableQuantity::Displacement, Variable::VariableType::vector, 2, NULL, {1,2}); + +#define _IFT_MagnetoElasticQuad_ll_Name "magnetoelasticquad_ll" +REGISTER_Element(MagnetoElasticQuad_ll) + + + + } // end namespace oofem diff --git a/src/mpm/termlibrary4.C b/src/mpm/termlibrary4.C index 5dce09c56..b705546f5 100644 --- a/src/mpm/termlibrary4.C +++ b/src/mpm/termlibrary4.C @@ -229,7 +229,7 @@ MagnetoElasticity_GradGrad_Term :: evaluate_lin (FloatMatrix& answer, MPElement& // Grad contains deformation gradient and magnetic field auto size = this->computeGradientFields(vGrad, vB, cell, gp->giveNaturalCoordinates(), gp->giveMaterialMode(), tstep); FloatMatrix D1,D2,D3,D4; - this->compute_lin_num(D1,D2,D3,D4, cell, gp, tstep); + // this->compute_lin_num(D1,D2,D3,D4, cell, gp, tstep); // auto cs = cell.giveCrossSection(); auto mcs = dynamic_cast (cs); diff --git a/src/oofemlib/internalstatetype.h b/src/oofemlib/internalstatetype.h index f3242fe53..68129b081 100644 --- a/src/oofemlib/internalstatetype.h +++ b/src/oofemlib/internalstatetype.h @@ -196,7 +196,10 @@ namespace oofem { ENUM_ITEM_WITH_VALUE(IST_VolumeFraction, 146) \ ENUM_ITEM_WITH_VALUE(IST_X_LCS, 147) /*Unit vector in local coordinate system in the x direction (usable for diagrams of internal forces for VTK export)*/ \ ENUM_ITEM_WITH_VALUE(IST_Y_LCS, 148) \ - ENUM_ITEM_WITH_VALUE(IST_Z_LCS, 149) + ENUM_ITEM_WITH_VALUE(IST_Z_LCS, 149) \ + ENUM_ITEM_WITH_VALUE(IST_MagneticPotential, 150) \ + ENUM_ITEM_WITH_VALUE(IST_MagneticFieldVector, 151) \ + ENUM_ITEM_WITH_VALUE(IST_MagneticInductionVector, 152) /** diff --git a/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp b/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp index 2f65a4517..ef4e8e770 100644 --- a/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp +++ b/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp @@ -2,150 +2,233 @@ namespace FTensor { - - // A(i,j,k,l) * B(m,n) double contraction + // A(i,j,k,l)*B(l) template - class Tensor4_times_Tensor1_double + int Dim3, char i, char j, char k, char l> + class Tensor4_times_Tensor1_single4 { Tensor4_Expr iterA; - Tensor1_Expr iterB; + Tensor1_Expr iterB; - public: - typename promote::V operator()(const int N1, const int N2, const int N3) const + template + typename promote::V + eval(const int N1, const int N2, const int N3, const Number &) const + { + return iterA(N1, N2, N3, Current_Dim - 1) * iterB(Current_Dim - 1) + + eval(N1, N2, N3, Number()); + } + typename promote::V + eval(const int N1, const int N2, const int N3, const Number<1> &) const { - typename promote::V result(0); - for(int xx = 0; xx < DimX; ++xx) - { - // Permutation is where the indices get checked. - result += Permutation4().eval( - iterA, N1, N2, N3, xx) - * iterB(xx); - } - return result; + return iterA(N1, N2, N3, 0) * iterB(0); } - Tensor4_times_Tensor1_double( - const Tensor4_Expr &iter_a, - const Tensor1_Expr &iter_b) - : iterA(iter_a), iterB(iter_b) + public: + Tensor4_times_Tensor1_single4( + const Tensor4_Expr &a, + const Tensor1_Expr &b) + : iterA(a), iterB(b) {} + typename promote::V operator()(const int N1, const int N2, const int N3) const + { + return eval(N1, N2, N3, Number()); + } + }; - - - // A(i,j,k,l)*B(l) - template - auto + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim2, i, j, k> operator*(const Tensor4_Expr &a, const Tensor1_Expr &b) { using TensorExpr - = Tensor4_times_Tensor1_double; - return Tensor3_Expr::V, Dim0, Dim1, Dim2, i, - j, k>(TensorExpr(a, b)); + = Tensor4_times_Tensor1_single4; + return Tensor3_Expr::V, Dim0, Dim1, Dim2, i, j, k>(TensorExpr(a, b)); } - + ////////////////////////////////////////////////////////////////////////////////////// // B(l) * A(i,j,k,l) - template - auto - operator*(const Tensor1_Expr &b, + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim2, i, j, k> + operator*(const Tensor1_Expr &b, const Tensor4_Expr &a) { - return a * b; + using TensorExpr + = Tensor4_times_Tensor1_single4; + return Tensor3_Expr::V, Dim0, Dim1, Dim2, i, j, k>(TensorExpr(a, b)); } - - - // A(i,j,k,l)*B(k) + ///////////////////////////////////////////////////////////////////////////////////// + // A(i,j,k,l)*B(k) template - auto + class Tensor4_times_Tensor1_single3 + { + Tensor4_Expr iterA; + Tensor1_Expr iterB; + + template + typename promote::V + eval(const int N1, const int N2, const int N3, const Number &) const + { + return iterA(N1, N2, Current_Dim - 1, N3) * iterB(Current_Dim - 1) + + eval(N1, N2, N3, Number()); + } + typename promote::V + eval(const int N1, const int N2, const int N3, const Number<1> &) const + { + return iterA(N1, N2, 0, N3) * iterB(0); + } + + public: + Tensor4_times_Tensor1_single3( + const Tensor4_Expr &a, + const Tensor1_Expr &b) + : iterA(a), iterB(b) + {} + typename promote::V operator()(const int N1, const int N2, const int N3) const + { + return eval(N1, N2, N3, Number()); + } + + }; + + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim3, i, j, l> operator*(const Tensor4_Expr &a, const Tensor1_Expr &b) { using TensorExpr - = Tensor4_times_Tensor1_double; - return Tensor3_Expr::V, Dim0, Dim1, Dim3, i, - j, l>(TensorExpr(a, b)); + = Tensor4_times_Tensor1_single3; + return Tensor3_Expr::V, Dim0, Dim1, Dim3, i, j, l>(TensorExpr(a, b)); } - - - // B(l) * A(i,j,k,l) - template - auto + // B(k) * A(i,j,k,l) + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim3, i, j, l> operator*(const Tensor1_Expr &b, - const Tensor4_Expr &a) + const Tensor4_Expr &a) { - return a * b; + using TensorExpr + = Tensor4_times_Tensor1_single3; + return Tensor3_Expr::V, Dim0, Dim1, Dim3, i, j, l>(TensorExpr(a, b)); } - - - // A(i,j,k,l)*B(j) + ////////////////////////////////////////////////////////////////////////////////////// + // A(i,j,k,l)*B(j) template - auto + class Tensor4_times_Tensor1_single2 + { + Tensor4_Expr iterA; + Tensor1_Expr iterB; + + template + typename promote::V + eval(const int N1, const int N2, const int N3, const Number &) const + { + return iterA(N1, Current_Dim - 1, N2, N3) * iterB(Current_Dim - 1) + + eval(N1, N2, N3, Number()); + } + typename promote::V + eval(const int N1, const int N2, const int N3, const Number<1> &) const + { + return iterA(N1, 0, N2, N3) * iterB(0); + } + + public: + Tensor4_times_Tensor1_single2( + const Tensor4_Expr &a, + const Tensor1_Expr &b) + : iterA(a), iterB(b) + {} + typename promote::V operator()(const int N1, const int N2, const int N3) const + { + return eval(N1, N2, N3, Number()); + } + + }; + + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim2, i, k, l> operator*(const Tensor4_Expr &a, const Tensor1_Expr &b) { using TensorExpr - = Tensor4_times_Tensor1_double; - return Tensor3_Expr::V, Dim0, Dim1, Dim3, i, - k, l>(TensorExpr(a, b)); + = Tensor4_times_Tensor1_single2; + return Tensor3_Expr::V, Dim0, Dim2, Dim3, i, k, l>(TensorExpr(a, b)); } - - - - // B(l) * A(i,j,k,l) - template - auto + // A(i,j,k,l)*B(j) + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim2, i, k, l> operator*(const Tensor1_Expr &b, - const Tensor4_Expr &a) + const Tensor4_Expr &a) { - return a * b; + using TensorExpr + = Tensor4_times_Tensor1_single2; + return Tensor3_Expr::V, Dim0, Dim2, Dim3, i, k, l>(TensorExpr(a, b)); } + ////////////////////////////////////////////////////////////////////////////////////// + // A(i,j,k,l)*B(i) + template + class Tensor4_times_Tensor1_single1 + { + Tensor4_Expr iterA; + Tensor1_Expr iterB; + template + typename promote::V + eval(const int N1, const int N2, const int N3, const Number &) const + { + return iterA(Current_Dim - 1, N1, N2, N3) * iterB(Current_Dim - 1) + + eval(N1, N2, N3, Number()); + } + typename promote::V + eval(const int N1, const int N2, const int N3, const Number<1> &) const + { + return iterA(0, N1, N2, N3) * iterB(0); + } + public: + Tensor4_times_Tensor1_single1( + const Tensor4_Expr &a, + const Tensor1_Expr &b) + : iterA(a), iterB(b) + {} + typename promote::V operator()(const int N1, const int N2, const int N3) const + { + return eval(N1, N2, N3, Number()); + } + + }; - - // A(i,j,k,l)*B(i) - template - auto + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim2, j, k, l> operator*(const Tensor4_Expr &a, const Tensor1_Expr &b) { using TensorExpr - = Tensor4_times_Tensor1_double; - return Tensor3_Expr::V, Dim1, Dim2, Dim3, j, - k, l>(TensorExpr(a, b)); + = Tensor4_times_Tensor1_single1; + return Tensor3_Expr::V, Dim0, Dim1, Dim3, j, k, l>(TensorExpr(a, b)); } - - - // B(i) * A(i,j,k,l) - template - auto + // *B(i) *A(i,j,k,l) + template + Tensor3_Expr, typename promote::V, Dim0, Dim1, Dim2, j, k, l> operator*(const Tensor1_Expr &b, - const Tensor4_Expr &a) + const Tensor4_Expr &a) { - return a * b; + using TensorExpr + = Tensor4_times_Tensor1_single1; + return Tensor3_Expr::V, Dim0, Dim1, Dim3, j, k, l>(TensorExpr(a, b)); } - + ////////////////////////////////////////////// + } diff --git a/src/oofemlib/unknowntype.h b/src/oofemlib/unknowntype.h index 4430d8562..a3175e7b6 100644 --- a/src/oofemlib/unknowntype.h +++ b/src/oofemlib/unknowntype.h @@ -51,7 +51,7 @@ namespace oofem { ENUM_ITEM_WITH_VALUE(DeplanationFunction, 16) \ ENUM_ITEM_WITH_VALUE(MacroSlipVector, 17) \ ENUM_ITEM_WITH_VALUE(ResidualForce, 18) \ - ENUM_ITEM_WITH_VALUE(magneticPotential, 19) + ENUM_ITEM_WITH_VALUE(MagneticPotential, 19) /** * Type representing particular unknown (its physical meaning). */ diff --git a/src/oofemlib/vtkbaseexportmodule.C b/src/oofemlib/vtkbaseexportmodule.C index f596005d0..b241202b6 100644 --- a/src/oofemlib/vtkbaseexportmodule.C +++ b/src/oofemlib/vtkbaseexportmodule.C @@ -493,7 +493,11 @@ VTKBaseExportModule::getNodalVariableFromPrimaryField(FloatArray &answer, DofMan answer.resize(3); } iState = IST_MacroSlipVector; - } else { + } else if ( type == MagneticPotential ) { + dofIDMask.followedBy(M_Pot); + iState = IST_MagneticPotential; + answer.resize(1); + } else { OOFEM_ERROR("unsupported unknownType %s", __UnknownTypeToString(type) ); } diff --git a/src/oofemlib/vtkxmlexportmodule.C b/src/oofemlib/vtkxmlexportmodule.C index 951160634..65617a25c 100644 --- a/src/oofemlib/vtkxmlexportmodule.C +++ b/src/oofemlib/vtkxmlexportmodule.C @@ -498,7 +498,7 @@ VTKXMLExportModule::giveDataHeaders(std::string &pointHeader, std::string &cellH if ( type == DisplacementVector || type == EigenVector || type == VelocityVector || type == DirectorField || type == MacroSlipVector || type == ResidualForce ) { vectors += __UnknownTypeToString(type); vectors.append(" "); - } else if ( type == FluxVector || type == PressureVector || type == Temperature || type == Humidity || type == DeplanationFunction ) { + } else if ( type == FluxVector || type == PressureVector || type == Temperature || type == Humidity || type == DeplanationFunction || type == MagneticPotential) { scalars += __UnknownTypeToString(type); scalars.append(" "); } else {