Expressions for Shell Element in TACS (#116)

* add new a2d shell code

* update the shell assemble frame

* add new mat concat expression

* update shell element a2d headers

* add test for the matrix rotate frame expression

* bug fixes and shell strain fixes

* a2d mat rotate frame tests pass!

* all shell internal deriv tests pass now

* working symMatRotateFrame

* fix forward comp

* reverse + forward ad compatibility

* updated a2d obj include

* CUDA compatibility in A2D

* update to use A2D_FUNCTION for cuda headers

* device compatible sqrt, etc other operations

* close #ifndef A2D_FUNCTION CUDA part

* clang format

* remove a2dmatconcat.h

* clang format v3

include/a2dcore.h

@@ -31,4 +31,11 @@
 #include "ad/a2dvecouter.h"
 #include "ad/a2dvecsum.h"
 
+// shell routines
+#include "ad/shell/a2dshellassembleframe.h"
+// #include "ad/shell/a2dmatconcat.h"
+#include "ad/shell/a2dmatrotateframe.h"
+#include "ad/shell/a2dshellstrain.h"
+#include "ad/shell/a2dsymmatrotateframe.h"
+
 #endif  //  A2D_CORE_H

include/a2ddefs.h

@@ -8,8 +8,14 @@
 template <typename T>
 using A2D_complex_t = std::complex<T>;
 
+// CUDA headers
 #ifndef A2D_FUNCTION
-#define A2D_FUNCTION  // A2D_FUNCTION does nothing in this scenario
+#ifdef __CUDACC__
+#define A2D_FUNCTION \
+  __host__ __device__  // A2D_FUNCTION does nothing in this scenario
+#else                  // not __CUDACC__
+#define A2D_FUNCTION
+#endif
 #endif
 
 namespace A2D {

include/ad/a2dgemm.h

@@ -119,7 +119,7 @@ class MatMatMultExpr {
         "Can't perform second order forward with first order objects");
     constexpr ADseed seed = conditional_value<ADseed, forder == ADorder::FIRST,
                                               ADseed::b, ADseed::p>::value;
-    if constexpr (adA == ADiffType::ACTIVE) {
+    if constexpr (adA == ADiffType::ACTIVE && adB == ADiffType::ACTIVE) {
       constexpr bool additive = true;
       MatMatMultCore<T, N, M, K, L, P, Q, opA, opB>(
           GetSeed<seed>::get_data(A), get_data(B), GetSeed<seed>::get_data(C));

include/ad/a2dgreenstrain.h

@@ -43,7 +43,9 @@ class MatGreenStrainExpr {
   static_assert(get_diff_order<Utype>::order == order,
                 "ADorder does not match");
 
-  A2D_FUNCTION MatGreenStrainExpr(Utype& Ux, Etype& E) : Ux(Ux), E(E) {}
+  A2D_FUNCTION MatGreenStrainExpr(Utype& Ux, Etype& E) : Ux(Ux), E(E) {
+    // printf("made mat green strain expression\n");
+  }
 
   A2D_FUNCTION void eval() {
     if constexpr (etype == GreenStrainType::LINEAR) {

include/ad/a2dmatsum.h

@@ -217,6 +217,8 @@ class MatSumScaleExpr {
     }
   }
 
+  A2D_FUNCTION void bzero() { B.bzero(); }
+
   A2D_FUNCTION void reverse() {
     constexpr ADseed seed = ADseed::b;
     if constexpr (adA == ADiffType::ACTIVE) {
@@ -237,6 +239,8 @@ class MatSumScaleExpr {
     }
   }
 
+  A2D_FUNCTION void hzero() { C.hzero(); }
+
   A2D_FUNCTION void hreverse() {
     constexpr ADseed seed = ADseed::h;
     if constexpr (adA == ADiffType::ACTIVE) {

include/ad/a2dmatvecmult.h

@@ -214,6 +214,131 @@ A2D_FUNCTION auto MatVecMult(const Atype& A, A2DObj<xtype>& x,
   return MatVecMultExpr<op, const Atype, A2DObj<xtype>, A2DObj<ytype>>(A, x, y);
 }
 
+// now define MatScale
+template <typename T, int M, int N>
+A2D_FUNCTION void MatScale(const T alpha, const Mat<T, M, N>& x,
+                           Mat<T, M, N>& y) {
+  MatScaleCore<T, M, N>(alpha, get_data(x), get_data(y));
+}
+
+template <class dtype, class Atype, class Btype>
+class MatScaleExpr {
+ public:
+  // Extract the numeric type to use
+  typedef typename get_object_numeric_type<dtype>::type T;
+
+  // Extract the dimensions of the underlying vectors
+  static constexpr int M = get_matrix_rows<Atype>::size;
+  static constexpr int N = get_matrix_columns<Atype>::size;
+  static constexpr int size = get_num_matrix_entries<Atype>::size;
+
+  // Get the differentiation order from the output
+  static constexpr ADorder order = get_diff_order<Btype>::order;
+
+  // Get the types of the matrices
+  static constexpr ADiffType add = get_diff_type<dtype>::diff_type;
+  static constexpr ADiffType adA = get_diff_type<Atype>::diff_type;
+
+  // Make sure the matrix dimensions are consistent
+  static_assert((get_a2d_object_type<Atype>::value ==
+                 get_a2d_object_type<Btype>::value),
+                "Matrices are not all of the same type");
+
+  A2D_FUNCTION MatScaleExpr(dtype alpha, Atype& A, Btype& B)
+      : alpha(alpha), A(A), B(B) {}
+
+  A2D_FUNCTION void eval() {
+    MatScaleCore<T, M, N>(get_data(alpha), get_data(A), get_data(B));
+  }
+
+  A2D_FUNCTION void bzero() { B.bzero(); }
+
+  template <ADorder forder>
+  A2D_FUNCTION void forward() {
+    constexpr ADseed seed = conditional_value<ADseed, forder == ADorder::FIRST,
+                                              ADseed::b, ADseed::p>::value;
+
+    if constexpr (add == ADiffType::ACTIVE && adA == ADiffType::ACTIVE) {
+      MatScaleCore<T, M, N>(GetSeed<seed>::get_data(alpha), get_data(A),
+                            GetSeed<seed>::get_data(B));
+      VecAddCore<T, size>(get_data(alpha), GetSeed<seed>::get_data(A),
+                          GetSeed<seed>::get_data(B));
+    } else if constexpr (add == ADiffType::ACTIVE) {
+      MatScaleCore<T, M, N>(GetSeed<seed>::get_data(alpha), get_data(A),
+                            GetSeed<seed>::get_data(B));
+    } else if constexpr (adA == ADiffType::ACTIVE) {
+      MatScaleCore<T, M, N>(get_data(alpha), GetSeed<seed>::get_data(A),
+                            GetSeed<seed>::get_data(B));
+    }
+  }
+  A2D_FUNCTION void reverse() {
+    constexpr ADseed seed = ADseed::b;
+    if constexpr (add == ADiffType::ACTIVE) {
+      GetSeed<seed>::get_data(alpha) +=
+          VecDotCore<T, size>(GetSeed<seed>::get_data(B), get_data(A));
+    }
+    if constexpr (adA == ADiffType::ACTIVE) {
+      VecAddCore<T, size>(get_data(alpha), GetSeed<seed>::get_data(B),
+                          GetSeed<seed>::get_data(A));
+    }
+  }
+
+  A2D_FUNCTION void hzero() { B.hzero(); }
+
+  A2D_FUNCTION void hreverse() {
+    if constexpr (add == ADiffType::ACTIVE) {
+      GetSeed<ADseed::h>::get_data(alpha) +=
+          VecDotCore<T, size>(GetSeed<ADseed::h>::get_data(B), get_data(A));
+    }
+    if constexpr (adA == ADiffType::ACTIVE) {
+      VecAddCore<T, size>(get_data(alpha), GetSeed<ADseed::h>::get_data(B),
+                          GetSeed<ADseed::h>::get_data(B));
+    }
+    if constexpr (add == ADiffType::ACTIVE && adA == ADiffType::ACTIVE) {
+      GetSeed<ADseed::h>::get_data(alpha) += VecDotCore<T, size>(
+          GetSeed<ADseed::b>::get_data(B), GetSeed<ADseed::p>::get_data(A));
+      VecAddCore<T, size>(GetSeed<ADseed::p>::get_data(alpha),
+                          GetSeed<ADseed::b>::get_data(B),
+                          GetSeed<ADseed::h>::get_data(A));
+    }
+  }
+
+  dtype alpha;
+  Atype& A;
+  Btype& B;
+};
+
+template <class T, class Atype, class Btype>
+A2D_FUNCTION auto MatScale(ADObj<T>& alpha, ADObj<Atype>& x, ADObj<Btype>& y) {
+  return MatScaleExpr<ADObj<T>&, ADObj<Atype>, ADObj<Atype>>(alpha, x, y);
+}
+
+template <class T, class Atype, class Btype>
+A2D_FUNCTION auto MatScale(const T alpha, ADObj<Atype>& x, ADObj<Btype>& y) {
+  return MatScaleExpr<const T, ADObj<Atype>, ADObj<Atype>>(alpha, x, y);
+}
+
+template <class T, class Atype, class Btype>
+A2D_FUNCTION auto MatScale(ADObj<T>& alpha, const Atype& x, ADObj<Btype>& y) {
+  return MatScaleExpr<ADObj<T>&, const Atype, ADObj<Atype>>(alpha, x, y);
+}
+
+template <class T, class Atype, class Btype>
+A2D_FUNCTION auto MatScale(A2DObj<T>& alpha, A2DObj<Atype>& x,
+                           A2DObj<Btype>& y) {
+  return MatScaleExpr<A2DObj<T>&, A2DObj<Atype>, A2DObj<Atype>>(alpha, x, y);
+}
+
+template <class T, class Atype, class Btype>
+A2D_FUNCTION auto MatScale(const T alpha, A2DObj<Atype>& x, A2DObj<Btype>& y) {
+  return MatScaleExpr<const T, A2DObj<Atype>, A2DObj<Atype>>(alpha, x, y);
+}
+
+template <class T, class Atype, class Btype>
+A2D_FUNCTION auto MatScale(A2DObj<T>& alpha, const Atype& x, A2DObj<Btype>& y) {
+  return MatScaleExpr<A2DObj<T>&, const Atype, A2DObj<Atype>>(alpha, x, y);
+}
+
 namespace Test {
 
 template <MatOp op, typename T, int N, int M, int K, int P>

include/ad/a2dobj.h

@@ -4,6 +4,7 @@
 #include "../a2ddefs.h"
 #include "a2dmat.h"
 #include "a2dvec.h"
+#include "adscalar.h"
 
 namespace A2D {
 
@@ -303,6 +304,11 @@ struct __get_matrix_rows<Mat<T, N, M>> {
   static constexpr int size = N;
 };
 
+template <template <typename, int> class SymMat, typename T, int N>
+struct __get_matrix_rows<SymMat<T, N>> {
+  static constexpr int size = N;
+};
+
 template <class T>
 struct get_matrix_rows : __get_matrix_rows<typename remove_a2dobj<T>::type> {
   static_assert(get_a2d_object_type<T>::value == ADObjType::MATRIX,
@@ -322,6 +328,11 @@ struct __get_matrix_columns<Mat<T, N, M>> {
   static constexpr int size = M;
 };
 
+template <template <typename, int> class SymMat, typename T, int N>
+struct __get_matrix_columns<SymMat<T, N>> {
+  static constexpr int size = N;
+};
+
 template <class T>
 struct get_matrix_columns
     : __get_matrix_columns<typename remove_a2dobj<T>::type> {
@@ -736,6 +747,155 @@ A2D_FUNCTION T* get_data(A2DObj<Vec<T, n>&>& vec) {
   return vec.value().get_data();
 }
 
+// new ADScalar get_data  (SPE)
+template <class T, int N>
+struct __is_numeric_type<ADScalar<T, N>> : std::is_floating_point<T> {};
+
+template <class T, int N>
+struct __is_numeric_type<ADScalar<std::complex<T>, N>>
+    : std::is_floating_point<T> {};
+
+template <int N>
+struct __get_object_numeric_type<ADScalar<double, N>> {
+  using type = ADScalar<double, N>;
+};
+
+template <int N>
+struct __get_object_numeric_type<ADScalar<std::complex<double>, N>> {
+  using type = ADScalar<std::complex<double>, N>;
+};
+
+template <typename T, int N,
+          std::enable_if_t<is_numeric_type<T>::value, bool> = true>
+ADScalar<T, N>& get_data(ADScalar<T, N>& value) {
+  return value;
+}
+
+template <typename T, int N,
+          std::enable_if_t<is_numeric_type<T>::value, bool> = true>
+const ADScalar<T, N>& get_data(const ADScalar<T, N>& value) {
+  return value;
+}
+
+template <typename T, int N,
+          std::enable_if_t<is_numeric_type<T>::value, bool> = true>
+ADScalar<T, N>& get_data(ADObj<ADScalar<T, N>>& value) {
+  return value.value();
+}
+
+template <typename T, int N,
+          std::enable_if_t<is_numeric_type<T>::value, bool> = true>
+const ADScalar<T, N>& get_data(const ADObj<ADScalar<T, N>>& value) {
+  return value.value();
+}
+
+template <typename T, int N,
+          std::enable_if_t<is_numeric_type<T>::value, bool> = true>
+ADScalar<T, N>& get_data(A2DObj<ADScalar<T, N>>& value) {
+  return value.value();
+}
+
+template <typename T, int N,
+          std::enable_if_t<is_numeric_type<T>::value, bool> = true>
+const ADScalar<T, N>& get_data(const A2DObj<ADScalar<T, N>>& value) {
+  return value.value();
+}
+
+/**
+ * @brief Get data pointers from objects
+ */
+template <typename T, int N, int m, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(Mat<ADScalar<T, N>, m, n>& mat) {
+  return mat.get_data();
+}
+
+template <typename T, int N, int m, int n>
+A2D_FUNCTION const ADScalar<T, N>* get_data(
+    const Mat<ADScalar<T, N>, m, n>& mat) {
+  return mat.get_data();
+}
+
+template <typename T, int N, int m, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(ADObj<Mat<ADScalar<T, N>, m, n>>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(A2DObj<Mat<ADScalar<T, N>, m, n>>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(ADObj<Mat<ADScalar<T, N>, m, n>&>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(A2DObj<Mat<ADScalar<T, N>, m, n>&>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m>
+A2D_FUNCTION ADScalar<T, N>* get_data(SymMat<ADScalar<T, N>, m>& mat) {
+  return mat.get_data();
+}
+
+template <typename T, int N, int m>
+A2D_FUNCTION const ADScalar<T, N>* get_data(
+    const SymMat<ADScalar<T, N>, m>& mat) {
+  return mat.get_data();
+}
+
+template <typename T, int N, int m>
+A2D_FUNCTION ADScalar<T, N>* get_data(ADObj<SymMat<ADScalar<T, N>, m>>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m>
+A2D_FUNCTION ADScalar<T, N>* get_data(A2DObj<SymMat<ADScalar<T, N>, m>>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m>
+A2D_FUNCTION ADScalar<T, N>* get_data(ADObj<SymMat<ADScalar<T, N>, m>&>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int m>
+A2D_FUNCTION ADScalar<T, N>* get_data(A2DObj<SymMat<ADScalar<T, N>, m>&>& mat) {
+  return mat.value().get_data();
+}
+
+template <typename T, int N, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(Vec<ADScalar<T, N>, n>& vec) {
+  return vec.get_data();
+}
+
+template <typename T, int N, int n>
+A2D_FUNCTION const ADScalar<T, N>* get_data(const Vec<ADScalar<T, N>, n>& vec) {
+  return vec.get_data();
+}
+
+template <typename T, int N, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(ADObj<Vec<ADScalar<T, N>, n>>& vec) {
+  return vec.value().get_data();
+}
+
+template <typename T, int N, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(A2DObj<Vec<ADScalar<T, N>, n>>& vec) {
+  return vec.value().get_data();
+}
+
+template <typename T, int N, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(ADObj<Vec<ADScalar<T, N>, n>&>& vec) {
+  return vec.value().get_data();
+}
+
+template <typename T, int N, int n>
+A2D_FUNCTION ADScalar<T, N>* get_data(A2DObj<Vec<ADScalar<T, N>, n>&>& vec) {
+  return vec.value().get_data();
+}
+
 }  // namespace A2D
 
 #endif  // A2D_OBJECTS_H