Improve elastodynamics perf 85 (#159)

* Add fixed arrays.

* Rewrite more material models.

* Clean up code.

* Fix a couple of bugs in Gucci model, clean up code.

* Add license and change include guard names.

Code/Source/svFSI/CMakeLists.txt

@@ -112,6 +112,7 @@ set(CSRCS
  ls.h ls.cpp
  main.cpp
  mat_fun.h mat_fun.cpp
+ mat_fun_carray.h mat_fun_carray.cpp
  mat_models.h mat_models.cpp
  mesh.h mesh.cpp
  nn.h nn.cpp

Code/Source/svFSI/Tensor4.h

@@ -34,7 +34,7 @@
 
 #include <cstring>
 
-#define Tensor4_check_enabled
+//#define Tensor4_check_enabled
 
 /// @brief The Tensor4 template class implements a simple interface to 4th order tensors.
 //

Code/Source/svFSI/main.cpp

@@ -257,6 +257,8 @@ void iterate_solution(Simulation* simulation)
  dmsg << "com_mod.eq[cEq].sym: " << com_mod.eq[cEq].sym;
  //simulation->com_mod.timer.set_time();
  #endif
+ //std::cout << "-------------------------------------" << std::endl;
+ //std::cout << "inner_count: " << inner_count << std::endl;
 
  auto istr = "_" + std::to_string(cTS) + "_" + std::to_string(inner_count);
  iEqOld = cEq;
@@ -328,6 +330,8 @@ void iterate_solution(Simulation* simulation)
  for (int iM = 0; iM < com_mod.nMsh; iM++) {
  eq_assem::global_eq_assem(com_mod, cep_mod, com_mod.msh[iM], Ag, Yg, Dg);
  }
+ com_mod.R.write("R_as"+ istr);
+ com_mod.Val.write("Val_as"+ istr);
 
  com_mod.Val.write("Val_as"+ istr);
  com_mod.R.write("R_as"+ istr);

Code/Source/svFSI/mat_fun.cpp

@@ -46,7 +46,6 @@
 namespace mat_fun {
 
 Array<int> t_ind;
-//static Array<int> t_ind;
 
 /// @brief Double dot product of 2 square matrices
 //
@@ -933,6 +932,37 @@ transpose(const Array<double>& A)
  return result;
 }
 
+void mat_mul6x3(const Array<double>& A, const Array<double>& B, Array<double>& C)
+{
+ #define mat_mul6x3_unroll 
+ #ifdef mat_mul6x3_unroll 
+ auto a = A.data();
+ auto b = B.data();
+ auto c = C.data();
+
+ c[0] = a[0]*b[0] + a[6]*b[1] + a[12]*b[2] + a[18]*b[3] + a[24]*b[4] + a[30]*b[5];
+ c[1] = a[1]*b[0] + a[7]*b[1] + a[13]*b[2] + a[19]*b[3] + a[25]*b[4] + a[31]*b[5];
+ c[2] = a[2]*b[0] + a[8]*b[1] + a[14]*b[2] + a[20]*b[3] + a[26]*b[4] + a[32]*b[5];
+ c[3] = a[3]*b[0] + a[9]*b[1] + a[15]*b[2] + a[21]*b[3] + a[27]*b[4] + a[33]*b[5];
+ c[4] = a[4]*b[0] + a[10]*b[1] + a[16]*b[2] + a[22]*b[3] + a[28]*b[4] + a[34]*b[5];
+ c[5] = a[5]*b[0] + a[11]*b[1] + a[17]*b[2] + a[23]*b[3] + a[29]*b[4] + a[35]*b[5];
+
+ c[6] = a[0]*b[6] + a[6]*b[7] + a[12]*b[8] + a[18]*b[9] + a[24]*b[10] + a[30]*b[11];
+ c[7] = a[1]*b[6] + a[7]*b[7] + a[13]*b[8] + a[19]*b[9] + a[25]*b[10] + a[31]*b[11];
+ c[8] = a[2]*b[6] + a[8]*b[7] + a[14]*b[8] + a[20]*b[9] + a[26]*b[10] + a[32]*b[11];
+ c[9] = a[3]*b[6] + a[9]*b[7] + a[15]*b[8] + a[21]*b[9] + a[27]*b[10] + a[33]*b[11];
+ c[10] = a[4]*b[6] + a[10]*b[7] + a[16]*b[8] + a[22]*b[9] + a[28]*b[10] + a[34]*b[11];
+ c[11] = a[5]*b[6] + a[11]*b[7] + a[17]*b[8] + a[23]*b[9] + a[29]*b[10] + a[35]*b[11];
+
+ c[12] = a[0]*b[12] + a[6]*b[13] + a[12]*b[14] + a[18]*b[15] + a[24]*b[16] + a[30]*b[17];
+ c[13] = a[1]*b[12] + a[7]*b[13] + a[13]*b[14] + a[19]*b[15] + a[25]*b[16] + a[31]*b[17];
+ c[14] = a[2]*b[12] + a[8]*b[13] + a[14]*b[14] + a[20]*b[15] + a[26]*b[16] + a[32]*b[17];
+ c[15] = a[3]*b[12] + a[9]*b[13] + a[15]*b[14] + a[21]*b[15] + a[27]*b[16] + a[33]*b[17];
+ c[16] = a[4]*b[12] + a[10]*b[13] + a[16]*b[14] + a[22]*b[15] + a[28]*b[16] + a[34]*b[17];
+ c[17] = a[5]*b[12] + a[11]*b[13] + a[17]*b[14] + a[23]*b[15] + a[29]*b[16] + a[35]*b[17];
+ #endif
+}
+
 };
 
 

Code/Source/svFSI/mat_fun.h

@@ -36,9 +36,6 @@
 #include "Tensor4.h"
 #include "Vector.h"
 
-#include "mat_fun_fixed.h"
-
-
 /// @brief The classes defined here duplicate the data structures in the 
 /// Fortran MATFUN module defined in MATFUN.f. 
 ///
@@ -61,6 +58,8 @@ namespace mat_fun {
  Vector<double> mat_mul(const Array<double>& A, const Vector<double>& v);
  Array<double> mat_mul(const Array<double>& A, const Array<double>& B);
  void mat_mul(const Array<double>& A, const Array<double>& B, Array<double>& result);
+ void mat_mul6x3(const Array<double>& A, const Array<double>& B, Array<double>& C);
+
  Array<double> mat_symm(const Array<double>& A, const int nd);
  Array<double> mat_symm_prod(const Vector<double>& u, const Vector<double>& v, const int nd);
  double mat_trace(const Array<double>& A, const int nd);

Code/Source/svFSI/mat_fun_carray.cpp

@@ -0,0 +1,68 @@
+/**
+ * Copyright (c) Stanford University, The Regents of the University of California, and others.
+ *
+ * All Rights Reserved.
+ *
+ * See Copyright-SimVascular.txt for additional details.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "mat_fun_carray.h"
+
+namespace mat_fun_carray {
+
+Array<int> t_ind;
+
+//----------
+// ten_init
+//----------
+// Initialize tensor index pointer
+//
+void ten_init(const int nd)
+{
+ if (t_ind.size() != 0) {
+ return;
+ }
+
+ int nn = pow(nd, 4);
+ t_ind.resize(4, nn);
+
+ int ii = 0;
+ for (int l = 0; l < nd; l++) {
+ for (int k = 0; k < nd; k++) {
+ for (int j = 0; j < nd; j++) {
+ for (int i = 0; i < nd; i++) {
+ t_ind(0,ii) = i;
+ t_ind(1,ii) = j;
+ t_ind(2,ii) = k;
+ t_ind(3,ii) = l;
+ ii = ii + 1;
+ }
+ }
+ }
+ }
+}
+
+};