*** empty log message ***

git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@362 305cdda6-5ce1-45b3-a98d-dfc68c8b3305

library/CMakeLists.txt

@@ -157,7 +157,6 @@ MACRO(ADD_NEKTAR_EXECUTABLE name sources)
 
     TARGET_LINK_LIBRARIES(${name}
         optimized LibUtilities debug LibUtilities-g
-        optimized StdRegions debug StdRegions-g
         optimized ${BOOST_THREAD_LIB} debug ${BOOST_THREAD_DEBUG_LIB})
         
     INSTALL(TARGETS ${name} RUNTIME DESTINATION ${NEKTAR_BIN_DIR} OPTIONAL

library/LibUtilities/CMakeLists.txt

@@ -35,7 +35,12 @@ SET(LibUtilityHeaders
 ./LinearAlgebra/NekMatrixImpl.hpp
 ./LinearAlgebra/NekMatrixMetadata.hpp
 ./LinearAlgebra/NekVector.hpp
+./LinearAlgebra/NekVectorFwd.hpp
+./LinearAlgebra/NekVectorConstantSized.hpp
+./LinearAlgebra/NekVectorVariableSized.hpp
 ./LinearAlgebra/NekVectorMetadata.hpp
+./LinearAlgebra/NekVectorTypeTraits.hpp
+./LinearAlgebra/NekVectorCommon.hpp
 ./LinearAlgebra/NekPoint.hpp
 ./ExpressionTemplates/Accumulator.hpp
 ./ExpressionTemplates/ArithmeticConcepts.hpp

library/LibUtilities/LinearAlgebra/NekDiagonalMatrix.hpp

@@ -79,7 +79,7 @@ namespace Nektar
                 std::copy(ptr, ptr+m_numberOfElements, begin());
             }
 
-            NekMatrix(unsigned int numberOfElements, DataType* ptr, MatrixDataHolderType t = eCopy) :
+            NekMatrix(unsigned int numberOfElements, DataType* ptr, PointerWrapper t = eCopy) :
                 m_numberOfElements(numberOfElements),
                 m_data()
             {
@@ -348,6 +348,9 @@ namespace Nektar
 
 /**
     $Log: NekDiagonalMatrix.hpp,v $
+    Revision 1.5  2007/01/29 01:30:19  bnelson
+    Removed memory manager requirements.
+
     Revision 1.4  2007/01/23 03:12:49  jfrazier
     Added more conditional compilation directives for expression templates.
 

library/LibUtilities/LinearAlgebra/NekFullMatrix.hpp

@@ -72,7 +72,7 @@ namespace Nektar
                 Initialize(rows, columns, DataType(0));
             }
 
-            NekMatrix(unsigned int rows, unsigned int columns, DataType* ptr, MatrixDataHolderType t = eCopy) :
+            NekMatrix(unsigned int rows, unsigned int columns, DataType* ptr, PointerWrapper t = eCopy) :
                 m_rows(rows),
                 m_columns(columns),
                 m_data()
@@ -163,7 +163,13 @@ namespace Nektar
         }
 #endif
 
-        ~NekMatrix() {}
+        ~NekMatrix() 
+        {
+            m_rows = 0;
+            m_columns = 0;
+            m_data.reset();
+        }
+
 
         unsigned int GetRows() const { return m_rows; }
         unsigned int GetColumns() const { return m_columns; }
@@ -353,6 +359,9 @@ namespace Nektar
 
 /**
     $Log: NekFullMatrix.hpp,v $
+    Revision 1.6  2007/02/04 04:27:43  bnelson
+    Updated linear systems to work with normal objects instead of only shared pointers.
+
     Revision 1.5  2007/01/29 01:30:20  bnelson
     Removed memory manager requirements.
 

library/LibUtilities/LinearAlgebra/NekLinAlgAlgorithms.hpp

@@ -73,7 +73,7 @@ namespace Nektar
         {
             for(unsigned int i = 0; i <= j-1; ++i)
             {
-                r(i,j) = q[i].dot(x[j]);
+                r(i,j) = q[i].Dot(x[j]);
             }
             
             VectorType y = x[j];
@@ -98,5 +98,8 @@ namespace Nektar
 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_LIN_ALG_ALGORITHMS_HPP
 
 /**
-    $Log: $
+    $Log: NekLinAlgAlgorithms.hpp,v $
+    Revision 1.1  2006/11/20 03:39:06  bnelson
+    Initial Revision
+
 **/

library/LibUtilities/LinearAlgebra/NekLinSys.hpp

@@ -54,17 +54,14 @@ namespace Nektar
         typedef NekVector<DataType, vectorDim, space> VectorType;
         typedef NekMatrix<DataType, eDiagonal, eNormal, space> MatrixType;
 
-        static VectorType Solve(const boost::shared_ptr<MatrixType>& A, const boost::shared_ptr<const VectorType>& b)
+        static void Solve(const boost::shared_ptr<MatrixType>& A, const VectorType& b, VectorType& x)
         {
-            ASSERTL0(A->GetColumns() == b->GetRows(), "ERROR: NekLinSys::Solve matrix columns must equal vector rows");
-            VectorType result(*b);
+            ASSERTL0(A->GetColumns() == b.GetRows(), "ERROR: NekLinSys::Solve matrix columns must equal vector rows");
 
             for(unsigned int i = 0; i < A->GetColumns(); ++i)
             {
-                result[i] = (*b)[i]/(*A)(i,i);
+                x[i] = b[i]/(*A)(i,i);
             }
-
-            return result;
         }
     };
 
@@ -74,53 +71,30 @@ namespace Nektar
         typedef NekMatrix<DataType, form, BlockType, space> MatrixType;
         typedef NekVector<DataType, vectorDim, space> VectorType;
 
-        static VectorType Solve(const boost::shared_ptr<MatrixType>& A, const boost::shared_ptr<const VectorType>& b)
+        static void Solve(const boost::shared_ptr<MatrixType>& A, const VectorType& b, VectorType& x)
         {
-            VectorType result(*b);
-            Lapack::dgetrs(A->GetRows(),A->GetColumns(),A->GetPtr().get(),result.GetPtr());
-            return result;
+            Lapack::dgetrs(A->GetRows(),A->GetColumns(),A->GetPtr().get(),x.GetPtr());
         }
     };
 
-    template<typename MatrixType, typename VectorType>
-    class LinearSystem
+    template<typename MatrixType>
+    class LinearSystem;
+
+    template<typename DataType, NekMatrixForm form, MatrixBlockType BlockType, unsigned int space>
+    class LinearSystem<NekMatrix<DataType, form, BlockType, space> >
     {
         public:
-            typedef LinearSystem<MatrixType, VectorType> ThisType;
-            typedef VectorType ResultType; 
+            typedef NekMatrix<DataType, form, BlockType, space> MatrixType;
+            typedef LinearSystem<MatrixType> ThisType;
 
         public:
-            LinearSystem(const boost::shared_ptr<MatrixType>& theA,
-                         const boost::shared_ptr<const VectorType>& theB) :
-                A(theA),
-                b(theB)
-            {
-            }
-
-            LinearSystem(const boost::shared_ptr<MatrixType>& theA,
-                         const VectorType& theB) :
-                A(theA),
-                b(&theB, DeleteNothing<VectorType>())
-            {
-            }
-
-            LinearSystem(MatrixType& theA,
-                         const boost::shared_ptr<const VectorType>& theB) :
-                A(&theA, DeleteNothing<MatrixType>()),
-                b(theB)
-            {
-            }
-
-            LinearSystem(MatrixType& theA,
-                         const VectorType& theB) :
-                A(&theA, DeleteNothing<MatrixType>()),
-                b(&theB, DeleteNothing<VectorType>())
+            explicit LinearSystem(const boost::shared_ptr<MatrixType>& theA) :
+                A(theA)
             {
             }
 
             LinearSystem(const ThisType& rhs) :
-                A(rhs.A),
-                b(rhs.b)
+                A(rhs.A)
             {
             }
 
@@ -134,20 +108,59 @@ namespace Nektar
             ~LinearSystem() {}
 
 
-            ResultType Solve() const
+            template<unsigned int dim>
+            NekVector<DataType, dim, space> Solve(boost::shared_ptr<NekVector<DataType, dim, space> > b) const
             {
-                return LinearSystemSolver<MatrixType, VectorType>::Solve(A, b);
+                typedef NekVector<DataType, dim, space> VectorType;
+                VectorType x(*b);
+                LinearSystemSolver<MatrixType, VectorType>::Solve(A, *b, x);
+                return x;
             }
 
+            template<unsigned int dim>
+            NekVector<DataType, dim, space> Solve(const NekVector<DataType, dim, space>& b) const
+            {
+                typedef NekVector<DataType, dim, space> VectorType;
+                VectorType x(b);
+                LinearSystemSolver<MatrixType, VectorType>::Solve(A, b, x);
+                return x;
+            }
+
+            //template<typename VectorType>
+            //VectorType Solve(const boost::shared_ptr<const VectorType>& b,
+            //                 const boost::shared_ptr<VectorType>& x) const
+            //{
+            //    return LinearSystemSolver<MatrixType, VectorType>::Solve(A, *b, *x);
+            //}
+
+            //template<typename VectorType>
+            //VectorType Solve(const boost::shared_ptr<const VectorType>& b,
+            //                 const VectorType& x) const
+            //{
+            //    return LinearSystemSolver<MatrixType, VectorType>::Solve(A, *b, x);
+            //}
+
+            //template<typename VectorType>
+            //VectorType Solve(const VectorType& b,
+            //                 const boost::shared_ptr<VectorType>& x) const
+            //{
+            //    return LinearSystemSolver<MatrixType, VectorType>::Solve(A, b, *x);
+            //}
+
+            //template<typename VectorType>
+            //VectorType Solve(const VectorType& b,
+            //                 const VectorType& x) const
+            //{
+            //    return LinearSystemSolver<MatrixType, VectorType>::Solve(A, b, *x);
+            //}
+
         private:
             void swap(ThisType& rhs)
             {
                 std::swap(A, rhs.A);
-                std::swap(b, rhs.b);
             }
 
             boost::shared_ptr<MatrixType> A;
-            boost::shared_ptr<const VectorType> b;
     };
 }
 

library/LibUtilities/LinearAlgebra/NekMatrixFwd.hpp

@@ -39,14 +39,13 @@
 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_MATRIX_FWD_HPP
 
 #include <LibUtilities/LinearAlgebra/MatrixBlockType.h>
+#include <LibUtilities/LinearAlgebra/PointerWrapper.h>
+
 #include <boost/shared_ptr.hpp>
 
 namespace Nektar
 {
-    /// \brief Specifies if the pointer passed to a NekMatrix is copied into an internal 
-    ///         representation or used directly.
-    enum MatrixDataHolderType { eWrapper, eCopy };
-    
+  
     /// \brief Enumeration used internally to determine if a pointer can be deleted or not.
     //enum MatrixDataDeltetableType { eDeletable, eNotDeletable };
     
@@ -59,7 +58,10 @@ namespace Nektar
 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_MATRIX_FWD_HPP
 
 /**
-    $Log: $
+    $Log: NekMatrixFwd.hpp,v $
+    Revision 1.5  2006/12/17 22:36:35  bnelson
+    Removed Macintosh line endings.
+
 **/
 
 

library/LibUtilities/LinearAlgebra/NekVectorCommon.hpp

+///////////////////////////////////////////////////////////////////////////////
+//
+// File: NekConstantSizedVector.hpp
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// 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.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: Constant sized and variable sized vectors have the same interface
+//              and performt he same algorithms - the only real difference is how
+//              the data is stored.  These methods allow use to use the same 
+//              algorithm regardless of data storage mechanism.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_VECTOR_COMMON_HPP
+#define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_VECTOR_COMMON_HPP
+
+#include <LibUtilities/LinearAlgebra/NekVectorFwd.hpp>
+#include <LibUtilities/LinearAlgebra/NekVectorTypeTraits.hpp>
+
+namespace Nektar
+{
+    template<typename DataType>
+    std::vector<DataType> fromString(const std::string& str)
+    {
+        std::vector<DataType> result;
+
+        try
+        {
+            typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
+            boost::char_separator<char> sep("(<,>) ");
+            tokenizer tokens(str, sep);
+            for( tokenizer::iterator strIter = tokens.begin(); strIter != tokens.end(); ++strIter)
+            {
+                result.push_back(boost::lexical_cast<DataType>(*strIter));
+            }
+        }
+
+        return result;
+    }
+
+    /// \todo Do the Norms with Blas where applicable.
+    template<typename DataType, unsigned int dim, unsigned int space>
+    DataType L1Norm(const NekVector<DataType, dim, space>& v)
+    {
+        typedef NekVector<DataType, dim, space> VectorType;
+
+        DataType result(0);
+        for(VectorType::const_iterator iter = v.begin(); iter != v.end(); ++iter)
+        {
+            result += NekVectorTypeTraits<DataType>::abs(*iter);
+        }
+
+        return result;
+    }
+    
+    template<typename DataType, unsigned int dim, unsigned int space>
+    DataType L2Norm(const NekVector<DataType, dim, space>& v)
+    {
+        typedef NekVector<DataType, dim, space> VectorType;
+
+        DataType result(0);
+        for(VectorType::const_iterator iter = v.begin(); iter != v.end(); ++iter)
+        {
+            DataType v = NekVectorTypeTraits<DataType>::abs(*iter);
+            result += v*v;
+        }
+        return sqrt(result);
+    }
+    
+    template<typename DataType, unsigned int dim, unsigned int space>
+    DataType InfinityNorm(const NekVector<DataType, dim, space>& v) 
+    {
+        DataType result = NekVectorTypeTraits<DataType>::abs(v[0]);
+        for(unsigned int i = 1; i < v.GetDimension(); ++i)
+        {
+            result = std::max(NekVectorTypeTraits<DataType>::abs(v[i]), result);
+        }
+        return result;
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    NekVector<DataType, dim, space> Negate(const NekVector<DataType, dim, space>& v) 
+    {
+        NekVector<DataType, dim, space> temp(v);
+        for(unsigned int i=0; i < temp.GetDimension(); ++i)
+        {
+            temp(i) = -temp(i);
+        }
+        return temp;
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    void PlusEqual(NekVector<DataType, dim, space>& lhs, const NekVector<DataType, dim, space>& rhs)
+    {
+        ASSERTL1(lhs.GetDimension() == rhs.GetDimension(), "Two vectors must have the same size in PlusEqual.")
+        
+        for(unsigned int i=0; i < lhs.GetDimension(); ++i)
+        {
+            lhs[i] += rhs[i];
+        }
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    void MinusEqual(NekVector<DataType, dim, space>& lhs, const NekVector<DataType, dim, space>& rhs)
+    {
+        ASSERTL1(lhs.GetDimension() == rhs.GetDimension(), "Two vectors must have the same size in MinusEqual.")
+        
+        for(unsigned int i=0; i < lhs.GetDimension(); ++i)
+        {
+            lhs[i] -= rhs[i];
+        }
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    void TimesEqual(NekVector<DataType, dim, space>& lhs, const DataType& rhs)
+    {
+        for(unsigned int i=0; i < lhs.GetDimension(); ++i)
+        {
+            lhs[i] *= rhs;
+        }
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    void DivideEqual(NekVector<DataType, dim, space>& lhs, const DataType& rhs)
+    {
+        for(unsigned int i=0; i < lhs.GetDimension(); ++i)
+        {
+            lhs[i] /= rhs;
+        }
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    DataType Magnitude(const NekVector<DataType, dim, space>& v) 
+    {
+        DataType result = DataType(0);
+
+        for(unsigned int i = 0; i < v.GetDimension(); ++i)
+        {
+            result += v[i]*v[i];
+        }
+        return sqrt(result);
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    DataType Dot(const NekVector<DataType, dim, space>& lhs, 
+                 const NekVector<DataType, dim, space>& rhs) 
+    {
+        ASSERTL1( lhs.GetDimension() == rhs.GetDimension(), "Dot, dimension of the two operands must be identical.");
+
+        DataType result = DataType(0);
+        for(unsigned int i = 0; i < lhs.GetDimension(); ++i)
+        {
+            result += lhs[i]*rhs[i];
+        }
+
+        return result;
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    void Normalize(NekVector<DataType, dim, space>& v)
+    {
+        DataType m = v.Magnitude();
+        if( m > DataType(0) )
+        {
+            v /= m;
+        }
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    NekVector<DataType, dim, space> Cross(const NekVector<DataType, dim, space>& lhs,
+                                          const NekVector<DataType, dim, space>& rhs)
+    {
+        BOOST_STATIC_ASSERT(dim==3 || dim == 0);
+        ASSERTL1(lhs.GetDimension() == 3 && rhs.GetDimension() == 3, "Cross is only valid for 3D vectors.");
+
+        DataType first = lhs.y()*rhs.z() - lhs.z()*rhs.y();
+        DataType second = lhs.z()*rhs.x() - lhs.x()*rhs.z();
+        DataType third = lhs.x()*rhs.y() - lhs.y()*rhs.x();
+
+        NekVector<DataType, dim, space> result(first, second, third);
+        return result;
+    }
+
+    template<typename DataType, unsigned int dim, unsigned int space>
+    std::string AsString(const NekVector<DataType, dim, space>& v)
+    {
+        unsigned int d = v.GetRows();
+        std::string result = "(";
+        for(unsigned int i = 0; i < d; ++i)
+        {
+            result += boost::lexical_cast<std::string>(v[i]);
+            if( i < dim-1 )
+            {
+                result += ", ";
+            }
+        }
+        result += ")";
+        return result;
+    }
+
+}
+
+#endif NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_VECTOR_COMMON_HPP
+

library/LibUtilities/LinearAlgebra/NekVectorConstantSized.hpp

+///////////////////////////////////////////////////////////////////////////////
+//
+// File: NekConstantSizedVector.hpp
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// 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.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: Generic N-Dimensional Vector.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_CONSTANT_SIZED_VECTOR_HPP
+#define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_NEK_CONSTANT_SIZED_VECTOR_HPP
+
+#include <LibUtilities/LinearAlgebra/NekVectorFwd.hpp>
+#include <LibUtilities/LinearAlgebra/NekVectorTypeTraits.hpp>
+#include <LibUtilities/LinearAlgebra/NekVectorCommon.hpp>
+
+#include <LibUtilities/ExpressionTemplates/ExpressionTemplates.hpp>
+#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
+
+
+#include <algorithm>
+#include <boost/static_assert.hpp>
+
+
+namespace Nektar
+{
+   
+    // \param DataType The type of data held by each element of the vector.
+    // \param dim The number of elements in the vector.  If set to 0, the vector
+    //            will have a variable number of elements.
+    // \param space The space of the vector.
+    template<typename DataType, unsigned int dim, unsigned int space>
+    class NekVector
+    {
+        public:
+            /// \brief Creates a constant sized vector with each element initialized to
+            ///        the default value for DataType.
+            NekVector();
+
+            /// \brief Creates a constant sized vector with each element set to a.
+            /// \param a The value to assign to each element of the vector.
+            NekVector(typename boost::call_traits<DataType>::const_reference a);
+
+            /// \brief Creates a vector from the elements in a delimited string.
+            /// \param vectorValues A string the the vector values.
+            /// 
+            /// 
+            explicit NekVector(const std::string& vectorValues);
+
+            NekVector(typename boost::call_traits<DataType>::const_reference x,
+                      typename boost::call_traits<DataType>::const_reference y,
+                      typename boost::call_traits<DataType>::const_reference z);
+
+            NekVector(typename boost::call_traits<DataType>::const_reference x,
+