1 // Copyright (C) 2022-2024 CEA, EDF
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 #include "InterpKernelDenseMatrix.hxx"
23 #include "InterpKernelException.hxx"
24 #include "VectorUtils.hxx"
28 namespace INTERP_KERNEL
31 DenseMatrixT<T>::DenseMatrixT() : nn(0), mm(0), v(nullptr) {}
34 DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
37 if (v) v[0] = nel>0 ? new T[nel] : nullptr;
38 for (i=1;i<n;i++) v[i] = v[i-1] + m;
42 DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m, const T &a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
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;
51 DenseMatrixT<T>::DenseMatrixT(mcIdType n, mcIdType m, const T *a) : nn(n), mm(m), v(n>0 ? new T*[n] : nullptr)
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++;
60 DenseMatrixT<T>::DenseMatrixT(const DenseMatrixT &rhs) : nn(rhs.nn), mm(rhs.mm), v(nn>0 ? new T*[nn] : nullptr)
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];
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
74 DenseMatrixT<T> & DenseMatrixT<T>::operator=(const DenseMatrixT<T> &rhs)
78 if (nn != rhs.nn || mm != rhs.mm) {
85 v = nn>0 ? new T*[nn] : nullptr;
87 if (v) v[0] = nel>0 ? new T[nel] : nullptr;
88 for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
90 for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
96 void DenseMatrixT<T>::resize(mcIdType newn, mcIdType newm)
99 if (newn != nn || newm != mm)
108 v = nn>0 ? new T*[nn] : nullptr;
110 if (v) v[0] = nel>0 ? new T[nel] : nullptr;
111 for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
116 void DenseMatrixT<T>::assign(mcIdType newn, mcIdType newm, const T& a)
119 if (newn != nn || newm != mm)
128 v = nn>0 ? new T*[nn] : nullptr;
130 if (v) v[0] = nel>0 ? new T[nel] : nullptr;
131 for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
133 for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = a;
137 DenseMatrixT<T>::~DenseMatrixT()
147 T Determinant22(const T *m)
148 { return m[0]*m[3]-m[1]*m[2]; }
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]);}
155 T DenseMatrixT<T>::determinant() const
160 return Determinant22(v[0]);
162 return Determinant33(v[0]);
163 THROW_IK_EXCEPTION("DenseMatrixT::determinant : only 1x1, 2x2 and 3x3 implemented !");
167 T DenseMatrixT<T>::toJacobian() const
170 return determinant();
171 const T *vPtr(this->v[0]);
175 return std::sqrt(vPtr[0]*vPtr[0] + vPtr[1]*vPtr[1]);
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]};
184 THROW_IK_EXCEPTION("DenseMatrixT::toJacobian : only 1x1, 2x1, 3x1, 3x2, 2x2 and 3x3 implemented !");