Adding Tensor4_times_Tensor1_single.hpp

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"

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 A, class B, class T, class U, int Dim0, int Dim1, int Dim2,
+            int Dim3, int Dim4, char i, char j, char k, char l,
+            char m, int DimA, int DimB, int DimC, int DimX, char a,
+            char b, char c, char x>
+  class Tensor4_times_Tensor1_double
+  {
+    Tensor4_Expr<A, T, Dim0, Dim1, Dim2, Dim3, i, j, k, l> iterA;
+    Tensor1_Expr<B, U, Dim4, m> iterB;
+
+  public:
+    typename promote<T, U>::V operator()(const int N1, const int N2, const int N3) const
+    {
+      typename promote<T, U>::V result(0);
+      for(int xx = 0; xx < DimX; ++xx)
+	{
+            // Permutation is where the indices get checked.
+            result += Permutation4<DimA, DimB, DimC, DimX, a, b, c, x>().eval(
+                        iterA, N1, N2, N3, xx)
+	              * iterB(xx);
+	}
+      return result;
+    }
+
+    Tensor4_times_Tensor1_double(
+      const Tensor4_Expr<A, T, Dim0, Dim1, Dim2, Dim3, i, j, k, l> &iter_a,
+      const Tensor1_Expr<B, U, Dim4, m> &iter_b)
+        : iterA(iter_a), iterB(iter_b)
+    {}
+  };
+
+
+
+  // A(i,j,k,l)*B(l)
+  template <class A, class B, class T, class U, int Dim0, int Dim1, int Dim2,
+            int Dim3, char i, char j, char k, char l>
+  auto
+  operator*(const Tensor4_Expr<A, T, Dim0, Dim1, Dim2, Dim3, i, j, k, l> &a,
+            const Tensor1_Expr<B, U, Dim3, l> &b)
+  {
+    using TensorExpr
+      = Tensor4_times_Tensor1_double<A, B, T, U, Dim0, Dim1, Dim2, Dim3, Dim3,
+				     i, j, k, l, l, Dim0, Dim1, Dim2,
+                                     Dim3, i, j, k, l>;
+    return Tensor3_Expr<TensorExpr, typename promote<T, U>::V, Dim0, Dim1, Dim2, i,
+                        j, k>(TensorExpr(a, b));
+  }
+
+   // A(i,j,k,l)*B(k)
+  template <class A, class B, class T, class U, int Dim0, int Dim1, int Dim2,
+            int Dim3, char i, char j, char k, char l>
+  auto
+  operator*(const Tensor4_Expr<A, T, Dim0, Dim1, Dim2, Dim3, i, j, k, l> &a,
+            const Tensor1_Expr<B, U, Dim2, k> &b)
+  {
+    using TensorExpr
+      = Tensor4_times_Tensor1_double<A, B, T, U, Dim0, Dim1, Dim2, Dim3, Dim2,
+				     i, j, k, l, k, Dim0, Dim1, Dim2,
+                                     Dim3, i, j, k, l>;
+    return Tensor3_Expr<TensorExpr, typename promote<T, U>::V, Dim0, Dim1, Dim3, i,
+                        j, l>(TensorExpr(a, b));
+  }
+
+   // A(i,j,k,l)*B(j)
+  template <class A, class B, class T, class U, int Dim0, int Dim1, int Dim2,
+            int Dim3, char i, char j, char k, char l>
+  auto
+  operator*(const Tensor4_Expr<A, T, Dim0, Dim1, Dim2, Dim3, i, j, k, l> &a,
+            const Tensor1_Expr<B, U, Dim1, j> &b)
+  {
+    using TensorExpr
+      = Tensor4_times_Tensor1_double<A, B, T, U, Dim0, Dim1, Dim2, Dim3, Dim1,
+				     i, j, k, l, j, Dim0, Dim1, Dim2,
+                                     Dim3, i, j, k, l>;
+    return Tensor3_Expr<TensorExpr, typename promote<T, U>::V, Dim0, Dim1, Dim3, i,
+                        k, l>(TensorExpr(a, b));
+  }
+
+
+   // A(i,j,k,l)*B(i)
+  template <class A, class B, class T, class U, int Dim0, int Dim1, int Dim2,
+            int Dim3, char i, char j, char k, char l>
+  auto
+  operator*(const Tensor4_Expr<A, T, Dim0, Dim1, Dim2, Dim3, i, j, k, l> &a,
+            const Tensor1_Expr<B, U, Dim0, i> &b)
+  {
+    using TensorExpr
+      = Tensor4_times_Tensor1_double<A, B, T, U, Dim0, Dim1, Dim2, Dim3, Dim0,
+				     i, j, k, l, i, Dim0, Dim1, Dim2,
+                                     Dim3, i, j, k, l>;
+    return Tensor3_Expr<TensorExpr, typename promote<T, U>::V, Dim1, Dim2, Dim3, j,
+                        k, l>(TensorExpr(a, b));
+  }
+
+
+
+}

src/oofemlib/tensor/tensor1.h

@@ -57,6 +57,13 @@ namespace oofem {
   {
   public:
     using Tensor1<double, 3> :: 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),
+      };
     }
 
   };

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<MagneticLoadElementInterface *>( 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;
+}//end namespace oofem;

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
+} // end namespace oofem

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;