From 7a7cf0e9c3cb352647b87f13696c475c64cef759 Mon Sep 17 00:00:00 2001 From: nitramkaroh Date: Thu, 22 Aug 2024 00:06:32 +0200 Subject: [PATCH] Adding Tensor4_times_Tensor1_single.hpp --- .../tensor/FTensor/Tensor4/Tensor4_Expr.hpp | 1 + .../Tensor4/Tensor4_times_Tensor1_single.hpp | 102 ++++++++++++++++++ src/oofemlib/tensor/tensor1.h | 13 ++- src/sm/Loads/hardmagneticboundarycondition.C | 32 +++--- src/sm/Loads/hardmagneticboundarycondition.h | 7 +- src/sm/Loads/magneticloadinterface.h | 5 +- 6 files changed, 139 insertions(+), 21 deletions(-) create mode 100644 src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp diff --git a/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_Expr.hpp b/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_Expr.hpp index 7cddc1575..e493ceb27 100644 --- a/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_Expr.hpp +++ b/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_Expr.hpp @@ -4,6 +4,7 @@ #include "Tensor4_minus_Tensor4.hpp" #include "Tensor4_plus_Tensor4.hpp" +#include "Tensor4_times_Tensor1_single.hpp" #include "Tensor4_times_Tensor2_single.hpp" #include "Tensor4_times_Tensor2_double.hpp" #include "Tensor4_times_Tensor2_symmetric.hpp" diff --git a/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp b/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp new file mode 100644 index 000000000..3ef41d6ea --- /dev/null +++ b/src/oofemlib/tensor/FTensor/Tensor4/Tensor4_times_Tensor1_single.hpp @@ -0,0 +1,102 @@ +#pragma once + +namespace FTensor +{ + + // A(i,j,k,l) * B(m,n) double contraction + template + class Tensor4_times_Tensor1_double + { + Tensor4_Expr iterA; + Tensor1_Expr iterB; + + public: + typename promote::V operator()(const int N1, const int N2, const int N3) 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; + } + + Tensor4_times_Tensor1_double( + const Tensor4_Expr &iter_a, + const Tensor1_Expr &iter_b) + : iterA(iter_a), iterB(iter_b) + {} + }; + + + + // A(i,j,k,l)*B(l) + template + auto + 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)); + } + + // A(i,j,k,l)*B(k) + template + auto + 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)); + } + + // A(i,j,k,l)*B(j) + template + auto + 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)); + } + + + // A(i,j,k,l)*B(i) + template + auto + 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)); + } + + + +} diff --git a/src/oofemlib/tensor/tensor1.h b/src/oofemlib/tensor/tensor1.h index e8d4a5568..bb290c111 100644 --- a/src/oofemlib/tensor/tensor1.h +++ b/src/oofemlib/tensor/tensor1.h @@ -57,6 +57,13 @@ namespace oofem { { public: using Tensor1 :: Tensor1; + + + Tensor1_3d() { + this->data [ 0 ] = 0; + this->data [ 1 ] = 0.; + this->data [ 2 ] = 0.; + } Tensor1_3d(const oofem::FloatArrayF<3> &array){ this->data[0] = array.at(1); this->data[1] = array.at(2); @@ -67,9 +74,9 @@ namespace oofem { { return { this->operator()(0), - this->operator()(1), - this->operator()(2), - }; + this->operator()(1), + this->operator()(2), + }; } }; diff --git a/src/sm/Loads/hardmagneticboundarycondition.C b/src/sm/Loads/hardmagneticboundarycondition.C index bfde3457c..b66791e73 100644 --- a/src/sm/Loads/hardmagneticboundarycondition.C +++ b/src/sm/Loads/hardmagneticboundarycondition.C @@ -160,12 +160,12 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition ); answer.clear(); + /* int order = 4; IntegrationRule *iRule = mfli->surfaceGiveIntegrationRule( order, iSurf ); - double a = 0; - + Tensor1_3d SigStarTimesCofN, Normal; - Tensor2_3d SigStar( load_level*load_level*sigma_star ); + Tensor2_3d SigStar( load_level*load_level*sigma_star ); Tensor2_3d F; FloatArray vF; @@ -189,11 +189,13 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition ); answer.plusProduct( N, SigStarTimesCofN.to_voigt_form(), gp->giveWeight() ); } + */ } void HardMagneticBoundaryCondition::computeTangentFromElement( FloatMatrix &answer, Element *e, int iSurf, TimeStep *tStep ) { + MagneticLoadElementInterface *mfli = static_cast( e->giveInterface( MagneticLoadElementInterfaceType ) ); if ( !mfli ) { @@ -215,20 +217,22 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition ); double nNodes = bNodes.giveSize(); FloatMatrix K; FloatMatrix testAnswer( 12, 12 ); + if ( e->giveSpatialDimension() == 3 ) { Tensor1_3d Normal; - Tensor2_3d SigStar( load_level * load_level * sigma_star ); + //@todo: this is not possible + Tensor2_3d SigStar( load_level * load_level * sigma_star ); Tensor2_3d F; Tensor3_3d SigStarFcrossN; FloatArray vF; - + for ( GaussPoint *gp : *iRule ) { FloatMatrix dN; // compute normal vector - FloatArray n, dxdk, dxde; - mfli->surfaceEvalNormalAt( n, dxdk, dxde, iSurf, gp, tStep ); + FloatArray n; + mfli->surfaceEvalNormalAt( n, iSurf, gp, tStep ); // compute surface N and B matrix FloatMatrix N, B; @@ -239,8 +243,8 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition ); mfli->surfaceEvalDeformationGradientAt( vF, iSurf, gp, tStep ); //compute the force vector - F = Tensor2_3d( vF ); - Normal = Tensor2_3d( n ); + F = Tensor2_3d( FloatArrayF< 9 >(vF) ); + Normal = Tensor1_3d( FloatArrayF< 3 >(n) ); auto [J, cofF] = F.compute_determinant_and_cofactor(); SigStarFcrossN( i_3, m_3, n_3 ) = SigStar( i_3, j_3) * F.compute_tensor_cross_product()( j_3, k_3, m_3, n_3) * Normal(k_3); @@ -256,18 +260,22 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition ); } else if ( e->giveSpatialDimension() == 2 ) { OOFEM_ERROR( "Magnetic boundary condition not implemented for 2D domains." ); } + } void HardMagneticBoundaryCondition::evaluateFreeSpaceStress() { + Tensor2_3d SigStar; - Tensor1_3d Bext( b_ext ); + FloatArrayF<3> bb( b_ext ); + Tensor1_3d Bext(bb); Tensor2_3d delta(1., 0., 0., 0., 1., 0., 0., 0., 1. ); - SigStar( i_3, j_3 ) = 1 / mu0 * (Bext(i_3)*Bext(j_3) - 0.5*Bext(p_3)*Bext(p_3)*delta(i_3,j_3)); + SigStar( i_3, j_3 ) = 1 / mu0 * (Bext(i_3) * Bext(j_3) - 0.5 * Bext(p_3) * Bext(p_3) * delta(i_3,j_3)); sigma_star = SigStar.to_voigt_form(); + } -}//end namespace oofem; \ No newline at end of file +}//end namespace oofem; diff --git a/src/sm/Loads/hardmagneticboundarycondition.h b/src/sm/Loads/hardmagneticboundarycondition.h index aaaba9f66..3f9ea400b 100644 --- a/src/sm/Loads/hardmagneticboundarycondition.h +++ b/src/sm/Loads/hardmagneticboundarycondition.h @@ -37,7 +37,8 @@ #include "activebc.h" #include "element.h" #include "inputrecord.h" - +#include "tensor/tensor1.h" +#include "tensor/tensor2.h" ///@name Input fields for HardMagneticBoundaryCondition //@{ #define _IFT_HardMagneticBoundaryCondition_Name "hardmagneticboundarycondition" @@ -61,7 +62,7 @@ class HardMagneticBoundaryCondition : public ActiveBoundaryCondition { double mu0; //vacuum permeability FloatArray b_ext; //external magnetic field in free space - FloatArray sigma_star; //precomputed free space stress + FloatArrayF<9> sigma_star; //precomputed free space stress int ltf_index; //index of load time function for applied load public: @@ -139,4 +140,4 @@ class HardMagneticBoundaryCondition : public ActiveBoundaryCondition { void evaluateFreeSpaceStress(); }; -} // end namespace oofem \ No newline at end of file +} // end namespace oofem diff --git a/src/sm/Loads/magneticloadinterface.h b/src/sm/Loads/magneticloadinterface.h index 54786d631..f60ff5f6a 100644 --- a/src/sm/Loads/magneticloadinterface.h +++ b/src/sm/Loads/magneticloadinterface.h @@ -60,15 +60,14 @@ class OOFEM_EXPORT MagneticLoadElementInterface : public Interface MagneticLoadElementInterface(Element *e); virtual ~MagneticLoadElementInterface(); - virtual const char *giveClassName() const { return "MagneticLoadElementInterface"; } + virtual const char *giveClassName() const override { return "MagneticLoadElementInterface"; } // private: virtual void surfaceEvalNmatrixAt( FloatMatrix &answer, int iSurf, GaussPoint *gp ) = 0; virtual void surfaceEvalBmatrixAt( FloatMatrix &answer, int iSurf, GaussPoint *gp ) = 0; - //virtual void surfaceEvaldNdxi( FloatMatrix &answer, int iSurf, GaussPoint *gp ) = 0; - virtual void surfaceEvalNormalAt(FloatArray &answer, FloatArray &dxdksi, FloatArray &dxdeta, int iSurf, GaussPoint *gp, TimeStep *tStep){;} + virtual void surfaceEvalNormalAt(FloatArray &answer, int iSurf, GaussPoint *gp, TimeStep *tStep){;} virtual void surfaceEvalDeformationGradientAt( FloatArray &answer, int isurf, GaussPoint *gp, TimeStep *tStep ) { ; } virtual IntegrationRule* surfaceGiveIntegrationRule(int order, int iSurf) = 0;