matrix::operator=

Apps/CMakeLists.txt

@@ -35,11 +35,11 @@ if(EXISTS ${PROJECT_SOURCE_DIR}/Elasticity)
   add_subdirectory(Elasticity/Linear)
   if(EXISTS ${PROJECT_SOURCE_DIR}/Elasticity/BeamEx)
     add_subdirectory(Elasticity/BeamEx)
-    add_subdirectory(Elasticity/BeamSim)
+    add_subdirectory(Elasticity/BeamEx/BeamSim)
   endif()
   if(EXISTS ${PROJECT_SOURCE_DIR}/Elasticity/FiniteDeformation)
     add_subdirectory(Elasticity/FiniteDeformation)
-    add_subdirectory(Elasticity/Nonlinear)
+    add_subdirectory(Elasticity/FiniteDeformation/Nonlinear)
   endif()
   if(EXISTS ${PROJECT_SOURCE_DIR}/ThermoElasticity)
     add_subdirectory(ThermoElasticity)

cmake/Modules/FindIFEM.cmake

@@ -27,14 +27,18 @@ if(NOT IFEM_AS_SUBMODULE)
     NAMES ${IFEM_BUILD_TYPE}/IFEM.h
     PATHS ${PROJECT_SOURCE_DIR}/..
           ${PROJECT_SOURCE_DIR}/../..
-          ${PROJECT_SOURCE_DIR}/../../.. NO_DEFAULT_PATHS)
+          ${PROJECT_SOURCE_DIR}/../../..
+          ${PROJECT_SOURCE_DIR}/../../../..
+          NO_DEFAULT_PATHS)
     set(IFEM_H_PATH ${IFEM_PATH}/${IFEM_BUILD_TYPE})
   if(NOT IFEM_PATH)
     find_path(IFEM_PATH
       NAMES IFEM.h
       PATHS ${PROJECT_SOURCE_DIR}/..
             ${PROJECT_SOURCE_DIR}/../..
-            ${PROJECT_SOURCE_DIR}/../../.. NO_DEFAULT_PATHS)
+            ${PROJECT_SOURCE_DIR}/../../..
+            ${PROJECT_SOURCE_DIR}/../../../..
+            NO_DEFAULT_PATHS)
     set(IFEM_H_PATH ${IFEM_PATH})
   endif(NOT IFEM_PATH)
 endif()
@@ -48,7 +52,8 @@ if(IFEM_AS_SUBMODULE)
   find_path(IFEM_PATH NAMES src/IFEM.h.in
                       PATHS ${PROJECT_SOURCE_DIR}/..
                             ${PROJECT_SOURCE_DIR}/../..
-                            ${PROJECT_SOURCE_DIR}/../../.. NO_DEFAULT_PATHS)
+                            ${PROJECT_SOURCE_DIR}/../../..
+                            ${PROJECT_SOURCE_DIR}/../../../.. NO_DEFAULT_PATHS)
   if(NOT IFEM_PATH)
     message(FATAL_ERROR "IFEM cannot be located, and we want it as a submodule")
   endif()

cmake/Scripts/IFEMTesting.cmake

@@ -41,7 +41,7 @@ macro(IFEM_add_test_app path workdir name)
 endmacro()
 
 macro(IFEM_add_unittests IFEM_PATH)
-  IFEM_add_test_app("${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C"
+  IFEM_add_test_app("${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C;${IFEM_PATH}/src/SIM/Test/*.C"
                     ${IFEM_PATH}
                     IFEM
                     ${IFEM_LIBRARIES} ${IFEM_DEPLIBS})

src/ASM/ASMmxBase.h

@@ -31,8 +31,7 @@ class ASMmxBase
 {
 protected:
   //! \brief The constructor sets the number of field variables.
-  //! \param[in] n_f1 Number of nodal variables in field 1
-  //! \param[in] n_f2 Number of nodal variables in field 2
+  //! \param[in] n_f Number of nodal variables in each field
   ASMmxBase(const std::vector<unsigned char>& n_f);
 
   //! \brief Initializes the patch level MADOF array.
@@ -48,24 +47,23 @@ class ASMmxBase
 			int basis = 0) const;
 
   //! \brief Injects nodal results for this patch into a global vector.
-  //! \param[in] globVec Global solution vector in DOF-order
-  //! \param[out] nodeVec Nodal result vector for this patch
-  //! \param[in] basis Which basis to extract the nodal values for
+  //! \param[out] globVec Global solution vector in DOF-order
+  //! \param[in] nodeVec Nodal result vector for this patch
+  //! \param[in] basis Which basis to inject the nodal values for
   void injectNodeVecMx(Vector& globVec, const Vector& nodeVec,
                        int basis = 0) const;
 
-  //! \brief Extract the primary solution field at the specified nodes.
+  //! \brief Extracts the primary solution field at the specified nodes.
   //! \param[out] sField Solution field
   //! \param[in] locSol Solution vector local to current patch
   //! \param[in] nodes 1-based local node numbers to extract solution for
   bool getSolutionMx(Matrix& sField, const Vector& locSol,
 		     const std::vector<int>& nodes) const;
 
 public:
-  static char geoBasis; //!< The basis representing the geometry
-
-  //! \brief Mixed formulation type
+  //! \brief Enum defining available mixed formulation types
   enum MixedType {
+    NONE = 0,                  //!< No fixed formulation
     FULL_CONT_RAISE_BASIS1,    //!< Full continuity, raise order and use as basis 1
     REDUCED_CONT_RAISE_BASIS1, //!< Reduced continuity, raise order and use as basis 1
     FULL_CONT_RAISE_BASIS2,    //!< Full continuity, raise order and use as basis 2
@@ -74,27 +72,29 @@ class ASMmxBase
   };
 
   static MixedType Type; //!< Type of mixed formulation used
+  static char  geoBasis; //!< 1-based index of basis representing the geometry
 
+protected:
   typedef std::vector<std::shared_ptr<Go::SplineSurface>> SurfaceVec; //!< Convenience type
   typedef std::vector<std::shared_ptr<Go::SplineVolume>> VolumeVec; //!< Convenience type
 
-  //! \brief Establish mixed bases
-  //! \param[in] surf The base basis to use.
-  //! \param[in] type The type of bases to establish.
-  //! \return Vector with bases.
+  //! \brief Establish mixed bases in 2D.
+  //! \param[in] surf The base basis to use
+  //! \param[in] type The type of bases to establish
+  //! \return Vector with bases
   static SurfaceVec establishBases(Go::SplineSurface* surf, MixedType type);
 
-  //! \brief Establish mixed bases
-  //! \param[in] vol The base basis to use.
-  //! \param[in] type The type of bases to establish.
-  //! \return Vector with bases.
+  //! \brief Establish mixed bases in 3D.
+  //! \param[in] vol The base basis to use
+  //! \param[in] type The type of bases to establish
+  //! \return Vector with bases
   static VolumeVec establishBases(Go::SplineVolume* vol, MixedType type);
 
 private:
   std::vector<int> MADOF; //!< Matrix of accumulated DOFs for this patch
 
 protected:
-  std::vector<size_t> nb;         //!< Number of basis functions in each basis
+  std::vector<size_t>        nb;  //!< Number of basis functions in each basis
   std::vector<unsigned char> nfx; //!< Number of fields on each basis
 };
 

src/ASM/IntegrandBase.h

@@ -59,7 +59,7 @@ class IntegrandBase : public Integrand
 
   //! \brief Defines the solution mode before the element assembly is started.
   virtual void setMode(SIM::SolutionMode mode) { m_mode = mode; }
-  //! \brief Obtain current solution mode
+  //! \brief Returns current solution mode.
   SIM::SolutionMode getMode() const { return m_mode; }
   //! \brief Initializes an integration parameter for the integrand.
   virtual void setIntegrationPrm(unsigned short int, double) {}
@@ -114,7 +114,7 @@ class IntegrandBase : public Integrand
                            const FiniteElement& fe,
                            const Vec3& X0, size_t nPt, LocalIntegral& elmInt);
   //! \brief Initializes current element for numerical integration (mixed).
-  //! \param[in] MNPC Nodal point correspondance for the bases
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] elem_sizes Size of each basis on the element
   //! \param[in] basis_sizes Size of each basis on the patch
   //! \param elmInt Local integral for element
@@ -129,7 +129,7 @@ class IntegrandBase : public Integrand
   virtual bool initElementBou(const std::vector<int>& MNPC,
                               LocalIntegral& elmInt);
   //! \brief Initializes current element for boundary integration (mixed).
-  //! \param[in] MNPC Nodal point correspondance for the bases
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
   //! \param[in] elem_sizes Size of each basis on the element
   //! \param[in] basis_sizes Size of each basis on the patch
   //! \param elmInt Local integral for element
@@ -153,7 +153,7 @@ class IntegrandBase : public Integrand
   //! \param[in] X Cartesian coordinates of current point
   //! \param[in] MNPC Nodal point correspondance for the basis function values
   virtual bool evalSol(Vector& s, const FiniteElement& fe, const Vec3& X,
-		       const std::vector<int>& MNPC) const;
+                       const std::vector<int>& MNPC) const;
 
   //! \brief Evaluates the secondary solution at a result point (mixed problem).
   //! \param[out] s The solution field values at current point
@@ -163,7 +163,7 @@ class IntegrandBase : public Integrand
   //! \param[in] elem_sizes Size of each basis on the element
   //! \param[in] basis_sizes Size of each basis on the patch
   virtual bool evalSol(Vector& s, const MxFiniteElement& fe, const Vec3& X,
-		       const std::vector<int>& MNPC,
+                       const std::vector<int>& MNPC,
                        const std::vector<size_t>& elem_sizes,
                        const std::vector<size_t>& basis_sizes) const;
 
@@ -172,21 +172,21 @@ class IntegrandBase : public Integrand
   //! \param[in] asol The analytical solution field (tensor field)
   //! \param[in] X Cartesian coordinates of current point
   virtual bool evalSol(Vector& s,
-		       const TensorFunc& asol, const Vec3& X) const;
+                       const TensorFunc& asol, const Vec3& X) const;
 
   //! \brief Evaluates the analytical secondary solution at a result point.
   //! \param[out] s The solution field values at current point
   //! \param[in] asol The analytical solution field (symmetric tensor field)
   //! \param[in] X Cartesian coordinates of current point
   virtual bool evalSol(Vector& s,
-		       const STensorFunc& asol, const Vec3& X) const;
+                       const STensorFunc& asol, const Vec3& X) const;
 
   //! \brief Evaluates the analytical secondary solution at a result point.
   //! \param[out] s The solution field values at current point
   //! \param[in] asol The analytical solution field (vector field)
   //! \param[in] X Cartesian coordinates of current point
   virtual bool evalSol(Vector& s,
-		       const VecFunc& asol, const Vec3& X) const;
+                       const VecFunc& asol, const Vec3& X) const;
 
   //! \brief Returns an evaluated principal direction vector field for plotting.
   virtual bool getPrincipalDir(Matrix&, size_t, size_t) const { return false; }
@@ -237,6 +237,8 @@ class IntegrandBase : public Integrand
 
   //! \brief Accesses the primary solution vector of current patch.
   Vector& getSolution(size_t n = 0) { return primsol[n]; }
+  //! \brief Accesses the primary solution vectors of current patch.
+  Vectors& getSolutions() { return primsol; }
 
   //! \brief Resets the primary solution vectors.
   void resetSolution();
@@ -258,7 +260,6 @@ class IntegrandBase : public Integrand
     return LinAlg::GENERAL_MATRIX;
   }
 
-protected:
   //! \brief Registers a vector to inject a named field into.
   //! \param[in] name Name of field
   //! \param[in] vec Vector to inject field into
@@ -300,6 +301,7 @@ class NormBase : public Integrand
   //! \brief Sets a vector of LocalIntegrals to be used during norm integration.
   void setLocalIntegrals(LintegralVec* elementNorms) { lints = elementNorms; }
 
+  using Integrand::getLocalIntegral;
   //! \brief Returns a local integral container for the element \a iEl.
   virtual LocalIntegral* getLocalIntegral(size_t, size_t iEl, bool) const;
 
@@ -401,6 +403,7 @@ class ForceBase : public Integrand
   //! \brief Initializes the integrand with the number of integration points.
   virtual void initIntegration(size_t, size_t) {}
 
+  using Integrand::getLocalIntegral;
   //! \brief Returns a local integral container for the element \a iEl.
   virtual LocalIntegral* getLocalIntegral(size_t, size_t iEl,
                                           bool = false) const;
@@ -415,7 +418,8 @@ class ForceBase : public Integrand
   { return false; }
 
   //! \brief Dummy implementation (only boundary integration is relevant).
-  virtual bool initElement(const std::vector<int>&, const std::vector<size_t>&,
+  virtual bool initElement(const std::vector<int>&,
+                           const std::vector<size_t>&,
                            const std::vector<size_t>&,
                            LocalIntegral&)
   { return false; }
@@ -424,7 +428,7 @@ class ForceBase : public Integrand
   virtual bool initElementBou(const std::vector<int>& MNPC,
                               LocalIntegral& elmInt);
   //! \brief Initializes current element for boundary integration (mixed).
-  virtual bool initElementBou(const std::vector<int>& MNPC1,
+  virtual bool initElementBou(const std::vector<int>& MNPC,
                               const std::vector<size_t>& elem_sizes,
                               const std::vector<size_t>& basis_sizes,
                               LocalIntegral& elmInt);

src/LinAlg/SAM.C

@@ -810,6 +810,8 @@ bool SAM::applyDirichlet (Vector& dofVec) const
       int ip = mpmceq[iceq-1];
       dofVec[idof] = ttcc[ip-1];
     }
+    else if (iceq == 0)
+      dofVec[idof] = 0.0;
   }
 
   return true;

src/LinAlg/Test/TestMatrix.C

@@ -17,25 +17,44 @@
 
 TEST(TestMatrix, AddBlock)
 {
-    utl::matrix<int> a(3,3);
-    int data_a[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
-    a.fill(data_a, 9);
-
-    utl::matrix<int> b(2,2);
-    int data_b[9] = {1, 2, 3, 4};
-    b.fill(data_b, 4);
-
-    a.addBlock(b, 2, 2, 2, false);
-
-    ASSERT_EQ(a(2,2), 7);
-    ASSERT_EQ(a(3,2), 10);
-    ASSERT_EQ(a(2,3), 14);
-    ASSERT_EQ(a(3,3), 17);
+  utl::matrix<int> a(3,3),b(2,2);
+  std::iota(a.begin(),a.end(),1);
+  std::iota(b.begin(),b.end(),1);
+
+  a.addBlock(b, 2, 2, 2, false);
+  ASSERT_EQ(a(2,2), 7);
+  ASSERT_EQ(a(3,2), 10);
+  ASSERT_EQ(a(2,3), 14);
+  ASSERT_EQ(a(3,3), 17);
+
+  a.addBlock(b, 1, 1, 1, true);
+  ASSERT_EQ(a(1,1), 2);
+  ASSERT_EQ(a(2,1), 5);
+  ASSERT_EQ(a(1,2), 6);
+  ASSERT_EQ(a(2,2), 11);
+}
 
-    a.addBlock(b, 1, 1, 1, true);
 
-    ASSERT_EQ(a(1,1), 2);
-    ASSERT_EQ(a(2,1), 5);
-    ASSERT_EQ(a(1,2), 6);
-    ASSERT_EQ(a(2,2), 11);
+TEST(TestMatrix, AddRows)
+{
+  utl::matrix<int> a(3,5);
+  std::iota(a.begin(),a.end(),1);
+  std::cout <<"A:"<< a;
+
+  a.expandCols(1);
+  std::cout <<"B:"<< a;
+  int fasit = 1;
+  for (size_t j = 1; j <= 5; j++)
+  {
+    for (size_t i = 1; i <= 3; i++, fasit++)
+      ASSERT_EQ(a(i,j), fasit);
+    ASSERT_EQ(a(4,j), 0);
+  }
+
+  a.expandCols(-2);
+  std::cout <<"C:"<< a;
+  fasit = 1;
+  for (size_t j = 1; j <= 5; j++, fasit++)
+    for (size_t i = 1; i <= 2; i++, fasit++)
+      ASSERT_EQ(a(i,j), fasit);
 }