src/INTERP_KERNEL/InterpKernelDenseMatrix.txx

   1 // Copyright (C) 2022  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #pragma once
  21
  22 #include "InterpKernelDenseMatrix.hxx"
  23 #include "InterpKernelException.hxx"
  24 #include "VectorUtils.hxx"
  25
  26 #include <cmath>
  27
  28 namespace INTERP_KERNEL
  29 {
  30   template <class T>
  31   DenseMatrixT<T>::DenseMatrixT() : nn(0), mm(0), v(nullptr) {}
  32
  33   template <class T>
  34   DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
  35   {
  36     mcIdType i,nel=m*n;
  37     if (v) v[0] = nel>0 ? new T[nel] : nullptr;
  38     for (i=1;i<n;i++) v[i] = v[i-1] + m;
  39   }
  40
  41   template <class T>
  42   DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m, const T &a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
  43   {
  44     mcIdType i,j,nel=m*n;
  45     if (v) v[0] = nel>0 ? new T[nel] : nullptr;
  46     for (i=1; i< n; i++) v[i] = v[i-1] + m;
  47     for (i=0; i< n; i++) for (j=0; j<m; j++) v[i][j] = a;
  48   }
  49
  50   template <class T>
  51   DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m, const T *a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
  52   {
  53     mcIdType i,j,nel=m*n;
  54     if (v) v[0] = nel>0 ? new T[nel] : nullptr;
  55     for (i=1; i< n; i++) v[i] = v[i-1] + m;
  56     for (i=0; i< n; i++) for (j=0; j<m; j++) v[i][j] = *a++;
  57   }
  58
  59   template <class T>
  60   DenseMatrixT<T>::DenseMatrixT(const DenseMatrixT &rhs) : nn(rhs.nn), mm(rhs.mm), v(nn>0 ? new T*[nn] : nullptr)
  61   {
  62     mcIdType i,j,nel=mm*nn;
  63     if (v) v[0] = nel>0 ? new T[nel] : nullptr;
  64     for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
  65     for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
  66   }
  67
  68   /*!
  69   * postcondition: normal assignment via copying has been performed
  70   * if matrix and rhs were different sizes, matrix
  71   * has been resized to match the size of rhs
  72   */
  73   template <class T>
  74   DenseMatrixT<T> & DenseMatrixT<T>::operator=(const DenseMatrixT<T> &rhs)
  75   {
  76     if (this != &rhs) {
  77       mcIdType i,j,nel;
  78       if (nn != rhs.nn || mm != rhs.mm) {
  79         if ( v ) {
  80           delete[] (v[0]);
  81           delete[] (v);
  82         }
  83         nn=rhs.nn;
  84         mm=rhs.mm;
  85         v = nn>0 ? new T*[nn] : nullptr;
  86         nel = mm*nn;
  87         if (v) v[0] = nel>0 ? new T[nel] : nullptr;
  88         for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
  89       }
  90       for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
  91     }
  92     return *this;
  93   }
  94
  95   template <class T>
  96   void DenseMatrixT<T>::resize(mcIdType newn, mcIdType newm)
  97   {
  98     mcIdType i,nel;
  99     if (newn != nn || newm != mm)
 100     {
 101       if ( v )
 102       {
 103         delete[] (v[0]);
 104         delete[] (v);
 105       }
 106       nn = newn;
 107       mm = newm;
 108       v = nn>0 ? new T*[nn] : nullptr;
 109       nel = mm*nn;
 110       if (v) v[0] = nel>0 ? new T[nel] : nullptr;
 111       for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
 112     }
 113   }
 114
 115   template <class T>
 116   void DenseMatrixT<T>::assign(mcIdType newn, mcIdType newm, const T& a)
 117   {
 118     mcIdType i,j,nel;
 119     if (newn != nn || newm != mm)
 120     {
 121       if ( v )
 122       {
 123         delete[] (v[0]);
 124         delete[] (v);
 125       }
 126       nn = newn;
 127       mm = newm;
 128       v = nn>0 ? new T*[nn] : nullptr;
 129       nel = mm*nn;
 130       if (v) v[0] = nel>0 ? new T[nel] : nullptr;
 131       for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
 132     }
 133     for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = a;
 134   }
 135
 136   template <class T>
 137   DenseMatrixT<T>::~DenseMatrixT()
 138   {
 139     if ( v )
 140     {
 141       delete[] (v[0]);
 142       delete[] (v);
 143     }
 144   }
 145
 146   template<class T>
 147   T Determinant22(const T *m)
 148   { return m[0]*m[3]-m[1]*m[2]; }
 149
 150   template<class T>
 151   T Determinant33(const T *m)
 152   { return m[0]*(m[4]*m[8]-m[7]*m[5])-m[1]*(m[3]*m[8]-m[6]*m[5])+m[2]*(m[3]*m[7]-m[6]*m[4]);}
 153
 154   template <class T>
 155   T DenseMatrixT<T>::determinant() const
 156   {
 157     if(nn==1 && mm==1)
 158       return v[0][0];
 159     if(nn==2 && mm==2)
 160       return Determinant22(v[0]);
 161     if(nn==3 && mm==3)
 162       return Determinant33(v[0]);
 163     THROW_IK_EXCEPTION("DenseMatrixT::determinant : only 1x1, 2x2 and 3x3 implemented !");
 164   }
 165
 166   template <class T>
 167   T DenseMatrixT<T>::toJacobian() const
 168   {
 169     if(nn == mm)
 170       return determinant();
 171     const T *vPtr(this->v[0]);
 172     if(nn==3 && mm==1)
 173       return norm(vPtr);
 174     if(nn==2 && mm==1)
 175       return std::sqrt(vPtr[0]*vPtr[0] + vPtr[1]*vPtr[1]);
 176     if(nn==3 && mm==2)
 177     {
 178       T tmp[3];
 179       T VA[3] = {v[0][0],v[1][0],v[2][0]};
 180       T VB[3] = {v[0][1],v[1][1],v[2][1]};
 181       cross(VA,VB,tmp);
 182       return norm(tmp);
 183     }
 184     THROW_IK_EXCEPTION("DenseMatrixT::toJacobian : only 1x1, 2x1, 3x1, 3x2, 2x2 and 3x3 implemented !");
 185   }
 186 }