diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index 41b160393..6c5ace82f 100644 --- a/src/ASM/ASMbase.C +++ b/src/ASM/ASMbase.C @@ -114,6 +114,8 @@ ASMbase::ASMbase (const ASMbase& patch) nLag = 0; // Lagrange multipliers are not copied myLMs.first = myLMs.second = 0; + + // why are these two added? Thought that by construction all vectors are empty neighbors.clear(); myLMTypes.clear(); } @@ -1039,19 +1041,29 @@ bool ASMbase::extractNodalVec (const Vector& globRes, Vector& nodeVec, nodeVec.reserve(nf*MLGN.size()); for (int inod : MLGN) { +#ifdef INDEX_CHECK + if (inod < 1 || (inod > ngnod && ngnod > 0)) + { + std::cerr <<" *** ASMbase::extractNodalVec: Global node "<< inod; + if (ngnod > 0) + std::cerr <<" is out of range [1,"<< ngnod <<"]."<< std::endl; + else + std::cerr <<" is out of range."<< std::endl; + return false; + } +#endif int idof = madof[inod-1] - 1; int jdof = madof[inod] - 1; + if (idof == jdof) + continue; // DOF-less node #ifdef INDEX_CHECK - bool ok = false; - if (inod < 1 || (inod > ngnod && ngnod > 0)) - std::cerr <<" *** ASMbase::extractNodalVec: Global node "<< inod - <<" is out of range [1,"<< std::abs(ngnod) <<"]."<< std::endl; - else if (jdof > (int)globRes.size()) - std::cerr <<" *** ASMbase::extractNodalVec: Global DOF "<< jdof - <<" is out of range [1,"<< globRes.size() <<"]."<< std::endl; - else - ok = true; - if (!ok) continue; + else if (idof < 0 || idof > jdof || jdof > (int)globRes.size()) + { + std::cerr <<" *** ASMbase::extractNodalVec: Global DOFs " + << idof+1 <<" "<< jdof + <<" out of range [1,"<< globRes.size() <<"]."<< std::endl; + return false; + } #endif nodeVec.insert(nodeVec.end(),globRes.ptr()+idof,globRes.ptr()+jdof); } @@ -1074,22 +1086,31 @@ bool ASMbase::injectNodalVec (const Vector& nodeVec, Vector& globVec, for (int inod : MLGN) if (basis == 0 || this->getNodeType(++i) == bType) { - int idof = madof[inod-1] - 1; - int ndof = madof[inod] - 1 - idof; #ifdef INDEX_CHECK - bool ok = false; if (inod < 1 || inod > (int)madof.size()) + { std::cerr <<" *** ASMbase::injectNodalVec: Node "<< inod <<" is out of range [1,"<< madof.size() <<"]."<< std::endl; - else if (ldof+ndof > nodeVec.size()) + return false; + } +#endif + int idof = madof[inod-1] - 1; + int ndof = madof[inod] - 1 - idof; + if (ndof == 0) + continue; // DOF-less node +#ifdef INDEX_CHECK + else if (ndof < 0 || ldof+ndof > nodeVec.size()) + { std::cerr <<" *** ASMbase::injectNodalVec: Local DOF "<< ldof+ndof <<" is out of range [1,"<< nodeVec.size() <<"]"<< std::endl; + return false; + } else if (idof+ndof > (int)globVec.size()) + { std::cerr <<" *** ASMbase::injectNodalVec: Global DOF "<< idof+ndof <<" is out of range [1,"<< globVec.size() <<"]"<< std::endl; - else - ok = true; - if (!ok) continue; + return false; + } #endif std::copy(nodeVec.begin()+ldof, nodeVec.begin()+ldof+ndof, globVec.begin()+idof); diff --git a/src/LinAlg/Test/TestMatrix.C b/src/LinAlg/Test/TestMatrix.C index a4ca8e134..b8abf566e 100644 --- a/src/LinAlg/Test/TestMatrix.C +++ b/src/LinAlg/Test/TestMatrix.C @@ -63,6 +63,35 @@ TEST(TestMatrix, AddRows) } +TEST(TestMatrix, Multiply) +{ + utl::vector u(14), v(9), x, y; + utl::matrix A(3,5); + + std::iota(u.begin(),u.end(),1.0); + std::iota(v.begin(),v.end(),1.0); + std::iota(A.begin(),A.end(),1.0); + + ASSERT_TRUE(A.multiply(u,x,1.0,0.0,false,3,4,1,2)); + ASSERT_TRUE(A.multiply(v,y,1.0,0.0,true,4,2)); + + EXPECT_FLOAT_EQ(x(3),370.0); + EXPECT_FLOAT_EQ(x(7),410.0); + EXPECT_FLOAT_EQ(x(11),450.0); + EXPECT_FLOAT_EQ(y(1),38.0); + EXPECT_FLOAT_EQ(y(3),83.0); + EXPECT_FLOAT_EQ(y(5),128.0); + EXPECT_FLOAT_EQ(y(7),173.0); + EXPECT_FLOAT_EQ(y(9),218.0); + + ASSERT_TRUE(A.multiply(u,x,1.0,-1.0,false,3,4,1,2)); + ASSERT_TRUE(A.multiply(v,y,1.0,-1.0,true,4,2)); + + EXPECT_FLOAT_EQ(x.sum(),0.0); + EXPECT_FLOAT_EQ(y.sum(),0.0); +} + + TEST(TestMatrix, Norm) { utl::matrix a(4,5); diff --git a/src/LinAlg/matrix.h b/src/LinAlg/matrix.h index 899b6f5af..010fbde75 100644 --- a/src/LinAlg/matrix.h +++ b/src/LinAlg/matrix.h @@ -92,7 +92,7 @@ namespace utl //! General utility classes and functions. } //! \brief Fill the vector with a scalar value. - void fill(const T& s) { std::fill(this->begin(),this->end(),s); } + void fill(T s) { std::fill(this->begin(),this->end(),s); } //! \brief Fill the vector with data from an array. void fill(const T* values, size_t n = 0) { @@ -125,7 +125,7 @@ namespace utl //! General utility classes and functions. //! \brief Subtract the given vector \b X from \a *this. vector& operator-=(const vector& X) { return this->add(X,T(-1)); } //! \brief Add the given vector \b X scaled by \a alfa to \a *this. - vector& add(const std::vector& X, T alfa = T(1)); + vector& add(const std::vector& X, const T& alfa = T(1)); //! \brief Perform \f${\bf Y} = \alpha{\bf Y} + (1-\alpha){\bf X} \f$ //! where \b Y = \a *this. @@ -317,7 +317,7 @@ namespace utl //! General utility classes and functions. void clear() { n[0] = n[1] = n[2] = 0; elem.clear(); } //! \brief Fill the matrix with a scalar value. - void fill(const T& s) { std::fill(elem.begin(),elem.end(),s); } + void fill(T s) { std::fill(elem.begin(),elem.end(),s); } //! \brief Fill the matrix with data from an array. void fill(const T* values, size_t n = 0) { elem.fill(values,n); } @@ -519,7 +519,7 @@ namespace utl //! General utility classes and functions. } //! \brief Add a scalar multiple of another matrix to a block of the matrix. - void addBlock(const matrix& block, const T& s, size_t r, size_t c, + void addBlock(const matrix& block, T s, size_t r, size_t c, bool transposed = false) { size_t nr = transposed ? block.cols() : block.rows(); @@ -533,7 +533,7 @@ namespace utl //! General utility classes and functions. } //! \brief Create a diagonal matrix. - matrix& diag(const T& d, size_t dim = 0) + matrix& diag(T d, size_t dim = 0) { if (dim > 0) this->resize(dim,dim,true); @@ -586,7 +586,7 @@ namespace utl //! General utility classes and functions. //! \brief Compute the inverse of a square matrix. //! \param[in] tol Division by zero tolerance //! \return Determinant of the matrix - T inverse(const T& tol = T(0)) + T inverse(T tol = T(0)) { T det = this->det(); if (det == T(-999)) @@ -626,7 +626,7 @@ namespace utl //! General utility classes and functions. //! \brief Check for symmetry. //! \param[in] tol Comparison tolerance - bool isSymmetric(const T& tol = T(0)) const + bool isSymmetric(T tol = T(0)) const { if (nrow != ncol) return false; @@ -673,7 +673,7 @@ namespace utl //! General utility classes and functions. */ matrix& multiply(const matrix& A, const matrix& B, bool transA = false, bool transB = false, - bool addTo = false, T alpha = T(1)); + bool addTo = false, const T& alpha = T(1)); /*! \brief Matrix-matrix multiplication. \details Performs one of the following operations (\b C = \a *this): @@ -723,8 +723,8 @@ namespace utl //! General utility classes and functions. -# \f$ {\bf Y} = {\alpha}{\bf A}^T {\bf X} + {\beta}{\bf Y}\f$ */ bool multiply(const std::vector& X, std::vector& Y, - T alpha, T beta, bool transA = false, - int stridex = 1, int stridey = 1, + const T& alpha, const T& beta = T(0), + bool transA = false, int stridex = 1, int stridey = 1, int ofsx = 0, int ofsy = 0) const; //! \brief Outer product between two vectors. @@ -1111,7 +1111,8 @@ namespace utl //! General utility classes and functions. } template<> inline - vector& vector::add(const std::vector& X, float alfa) + vector& vector::add(const std::vector& X, + const float& alfa) { size_t n = this->size() < X.size() ? this->size() : X.size(); cblas_saxpy(n,alfa,&X.front(),1,this->ptr(),1); @@ -1119,7 +1120,8 @@ namespace utl //! General utility classes and functions. } template<> inline - vector& vector::add(const std::vector& X, double alfa) + vector& vector::add(const std::vector& X, + const double& alfa) { size_t n = this->size() < X.size() ? this->size() : X.size(); cblas_daxpy(n,alfa,&X.front(),1,this->ptr(),1); @@ -1197,11 +1199,18 @@ namespace utl //! General utility classes and functions. template<> inline bool matrix::multiply(const std::vector& X, std::vector& Y, - float alpha, float beta, bool transA, - int stridex, int stridey, int ofsx, - int ofsy) const + const float& alpha, const float& beta, + bool transA, int stridex, int stridey, + int ofsx, int ofsy) const { - if (!this->compatible(X,transA)) return false; + if (ofsx == 0 && stridex == 1) + if (!this->compatible(X,transA)) return false; + + if (beta == 0.0f) + { + Y.resize(ofsy + 1 + ((transA ? ncol : nrow)-1)*stridey); + if (stridey > 0) std::fill(Y.begin(),Y.end(),0.0f); + } cblas_sgemv(CblasColMajor, transA ? CblasTrans : CblasNoTrans, @@ -1216,11 +1225,18 @@ namespace utl //! General utility classes and functions. template<> inline bool matrix::multiply(const std::vector& X, std::vector& Y, - double alpha, double beta, bool transA, - int stridex, int stridey, + const double& alpha, const double& beta, + bool transA, int stridex, int stridey, int ofsx, int ofsy) const { - if (!this->compatible(X,transA)) return false; + if (ofsx == 0 && stridex == 1) + if (!this->compatible(X,transA)) return false; + + if (beta == 0.0) + { + Y.resize(ofsy + 1 + ((transA ? ncol : nrow)-1)*stridey); + if (stridey > 0) std::fill(Y.begin(),Y.end(),0.0); + } cblas_dgemv(CblasColMajor, transA ? CblasTrans : CblasNoTrans, @@ -1236,20 +1252,23 @@ namespace utl //! General utility classes and functions. matrix& matrix::multiply(const matrix& A, const matrix& B, bool transA, bool transB, - bool addTo, float alpha) + bool addTo, const float& alpha) { size_t M, N, K; - if (!this->compatible(A,B,transA,transB,M,N,K)) return *this; - if (!addTo) this->resize(M,N); - - cblas_sgemm(CblasColMajor, - transA ? CblasTrans : CblasNoTrans, - transB ? CblasTrans : CblasNoTrans, - M, N, K, alpha, - A.ptr(), A.nrow, - B.ptr(), B.nrow, - addTo ? 1.0f : 0.0f, - this->ptr(), nrow); + if (!this->compatible(A,B,transA,transB,M,N,K)) + this->clear(); + else if (!addTo) + this->resize(M,N); + + if (!this->empty()) + cblas_sgemm(CblasColMajor, + transA ? CblasTrans : CblasNoTrans, + transB ? CblasTrans : CblasNoTrans, + M, N, K, alpha, + A.ptr(), A.nrow, + B.ptr(), B.nrow, + addTo ? 1.0f : 0.0f, + this->ptr(), nrow); return *this; } @@ -1258,20 +1277,23 @@ namespace utl //! General utility classes and functions. matrix& matrix::multiply(const matrix& A, const matrix& B, bool transA, bool transB, - bool addTo, double alpha) + bool addTo, const double& alpha) { size_t M, N, K; - if (!this->compatible(A,B,transA,transB,M,N,K)) return *this; - if (!addTo) this->resize(M,N); - - cblas_dgemm(CblasColMajor, - transA ? CblasTrans : CblasNoTrans, - transB ? CblasTrans : CblasNoTrans, - M, N, K, alpha, - A.ptr(), A.nrow, - B.ptr(), B.nrow, - addTo ? 1.0 : 0.0, - this->ptr(), nrow); + if (!this->compatible(A,B,transA,transB,M,N,K)) + this->clear(); + else if (!addTo) + this->resize(M,N); + + if (!this->empty()) + cblas_dgemm(CblasColMajor, + transA ? CblasTrans : CblasNoTrans, + transB ? CblasTrans : CblasNoTrans, + M, N, K, alpha, + A.ptr(), A.nrow, + B.ptr(), B.nrow, + addTo ? 1.0 : 0.0, + this->ptr(), nrow); return *this; } @@ -1416,7 +1438,7 @@ namespace utl //! General utility classes and functions. template<> inline bool matrix3d::multiplyMat(const std::vector& A, - const matrix3d& B, bool addTo) + const matrix3d& B, bool addTo) { size_t M, N, K; if (!this->compatible(A,B,M,N,K)) return false; @@ -1546,7 +1568,7 @@ namespace utl //! General utility classes and functions. } template inline - vector& vector::add(const std::vector& X, T alfa) + vector& vector::add(const std::vector& X, const T& alfa) { T* p = this->ptr(); const T* q = &X.front(); @@ -1556,7 +1578,7 @@ namespace utl //! General utility classes and functions. } template inline - matrixBase& matrixBase::add(const matrixBase& A, T alfa) + matrixBase& matrixBase::add(const matrixBase& A, const T& alfa) { T* p = this->ptr(); const T* q = A.ptr(); @@ -1596,20 +1618,28 @@ namespace utl //! General utility classes and functions. template inline bool matrix::multiply(const std::vector& X, std::vector& Y, - T alpha, T beta, bool transA, - int stridex, int stridey, + const T& alpha, const T& beta, + bool transA, int stridex, int stridey, int ofsx, int ofsy) const { - if (!this->compatible(X,transA)) return false; + if (ofsx == 0 && stridex == 1) + if (!this->compatible(X,transA)) return false; - for (size_t i = ofsy; i < Y.size(); i += stridey) { + if (beta == T(0)) + { + Y.resize(ofsy + 1 + ((transA ? ncol : nrow)-1)*stridey); + std::fill(Y.begin(),Y.end(),T(0)); + } + else for (size_t i = ofsy; i < Y.size(); i += stridey) Y[i] *= beta; - for (size_t j = ofsx; j < X.size(); j += stridex) + + size_t a, b, i, j; + for (a = 1, i = ofsy; i < Y.size(); a++, i += stridey) + for (b = 1, j = ofsx; j < X.size(); b++, j += stridex) if (transA) - Y[i] += alpha * THIS(j+1,i+1) * X[j]; + Y[i] += alpha * THIS(b,a) * X[j]; else - Y[i] += alpha * THIS(i+1,j+1) * X[j]; - } + Y[i] += alpha * THIS(a,b) * X[j]; return true; } @@ -1618,11 +1648,16 @@ namespace utl //! General utility classes and functions. matrix& matrix::multiply(const matrix& A, const matrix& B, bool transA, bool transB, bool addTo, - T alpha) + const T& alpha) { size_t M, N, K; - if (!this->compatible(A,B,transA,transB,M,N,K)) return *this; - if (!addTo) this->resize(M,N,true); + if (!this->compatible(A,B,transA,transB,M,N,K)) + { + this->clear(); + M = 0; + } + else if (!addTo) + this->resize(M,N,true); for (size_t i = 1; i <= M; i++) for (size_t j = 1; j <= N; j++) @@ -1749,7 +1784,7 @@ namespace utl //! General utility classes and functions. //! matrices when they contain terms that are numerically zero, except for //! some round-off noise. The value of the global variable \a zero_print_tol //! is used as a tolerance in this method. - template inline T trunc(const T& v) + template inline T trunc(T v) { return v > T(zero_print_tol) || v < T(-zero_print_tol) || std::isnan(v) ? v : T(0); diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index a471636d5..062b50b72 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -1623,11 +1623,15 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, ssol.clear(); - size_t i, j, n; - size_t ngNodes = this->getNoNodes(1); - - Matrix values; + size_t ngNodes = this->fieldProjections() ? 0 : this->getNoNodes(1); Vector count(myModel.size() > 1 ? ngNodes : 0); + Matrix values; + + if (this->fieldProjections()) + // No nodal averaging - full (potentially discontinuous) representation + for (ASMbase* pch : myModel) + if (!pch->empty()) + ngNodes += pch->getNoProjectionNodes(); if (pMethod == SIMoptions::DGL2 || pMethod == SIMoptions::CGL2) const_cast(this)->setQuadratureRule(opt.nGauss[1]); @@ -1639,17 +1643,7 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, myProblem->initIntegration(time,psol); } - // no nodal averaging - full (potentially discontinuous) representation - if (this->fieldProjections()) { - ngNodes = 0; - for (i = 0; i < myModel.size(); i++) - { - if (myModel[i]->empty()) continue; // skip empty patches - ngNodes += myModel[i]->getNoProjectionNodes(); - } - } - - size_t ofs = 0; + size_t i, ofs = 0; for (i = 0; i < myModel.size(); i++) { if (myModel[i]->empty()) continue; // skip empty patches @@ -1709,7 +1703,7 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, case SIMoptions::LEASTSQ: if (msgLevel > 1 && i == 0) - IFEM::cout <<"\tLeast squares projection"<< std::endl; + IFEM::cout <<"\tLeast squares projection"<< std::endl; ok = myModel[i]->evalSolution(values,*myProblem,nullptr,'W'); break; @@ -1741,13 +1735,13 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, // Nodal averaging for nodes referred to by two or more patches // (these are typically the interface nodes) - for (n = 1; n <= nNodes; n++) + for (size_t n = 1; n <= nNodes; n++) if (count.empty()) ssol.fillColumn(myModel[i]->getNodeID(n),values.getColumn(n)); else { int inod = myModel[i]->getNodeID(n); - for (j = 1; j <= nComps; j++) + for (size_t j = 1; j <= nComps; j++) ssol(j,inod) += values(j,n); count(inod) ++; } @@ -1755,11 +1749,10 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, } // Divide through by count(n) to get the nodal average at the interface nodes - if (!this->fieldProjections()) - for (n = 1; n <= count.size(); n++) - if (count(n) > 1.0) - for (j = 1; j <= ssol.rows(); j++) - ssol(j,n) /= count(n); + for (size_t n = 1; n <= count.size(); n++) + if (count(n) > 1.0) + for (size_t j = 1; j <= ssol.rows(); j++) + ssol(j,n) /= count(n); return true; } @@ -1796,7 +1789,8 @@ bool SIMbase::project (Vector& values, const FunctionBase* f, } if (nFields <= (int)f->dim()) - ok &= myModel[j]->injectNodeVec(loc_values,values,f->dim(),basis); + ok &= this->injectPatchSolution(values,loc_values, + myModel[j],f->dim(),basis); else if (f->dim() > 1) { std::cerr <<" *** SIMbase::project: Cannot interleave non-scalar function" @@ -1892,8 +1886,7 @@ size_t SIMbase::extractPatchSolution (const Vector& sol, Vector& vec, { if (!pch || sol.empty()) return 0; - if (basis != 0 && nndof != 0 && - nndof != this->getNoFields(basis) && this->getNoFields(2) > 0) + if (basis && nndof != pch->getNoFields(basis) && pch->getNoFields(2) > 0) { // Need an additional MADOF const IntVec& madof = this->getMADOF(basis,nndof); @@ -1912,12 +1905,11 @@ bool SIMbase::injectPatchSolution (Vector& sol, const Vector& vec, { if (!pch) return false; - if (basis > 0 && nndof > 0 && - nndof != this->getNoFields(basis) && this->getNoFields(2) > 0) + if (basis && nndof != pch->getNoFields(basis) && pch->getNoFields(2) > 0) // Need an additional MADOF return pch->injectNodalVec(vec,sol,this->getMADOF(basis,nndof),basis); else - return pch->injectNodeVec(vec,sol,nndof); + return pch->injectNodeVec(vec,sol,nndof,basis); } @@ -1944,9 +1936,9 @@ bool SIMbase::extractPatchElmRes (const Matrix& glbRes, Matrix& elRes, bool SIMbase::setPatchMaterial (size_t patch) { - for (PropertyVec::const_iterator p = myProps.begin(); p != myProps.end(); ++p) - if (p->pcode == Property::MATERIAL && p->patch == patch) - return this->initMaterial(p->pindx); + for (const Property& p : myProps) + if (p.pcode == Property::MATERIAL && p.patch == patch) + return this->initMaterial(p.pindx); return false; } @@ -1962,22 +1954,36 @@ bool SIMbase::addMADOF (unsigned char basis, unsigned char nndof, bool other) madof.resize(this->getNoNodes()+1,0); if (madof.size() < 2) return false; + int ierr = 0; char nType = basis <= 1 ? 'D' : 'P' + basis-2; - for (size_t i = 0; i < myModel.size(); i++) - for (size_t j = 0; j < myModel[i]->getNoNodes(); j++) + for (const ASMbase* pch : myModel) + for (size_t inod = 1; inod <= pch->getNoNodes(); inod++) { - int n = myModel[i]->getNodeID(j+1); - if (n > 0 && myModel[i]->getNodeType(j+1) == nType) - madof[n] = nndof; - else if (n > 0 && other) - madof[n] = myModel[i]->getNodalDOFs(j+1); + int node = pch->getNodeID(inod); + if (node > 0) + { + int nndofs = 0; + if (pch->getNodeType(inod) == nType) + nndofs = nndof; + else if (other) + nndofs = pch->getNodalDOFs(inod); + if (nndofs > 0 && madof[node] == 0) + madof[node] = nndofs; + else if (madof[node] != nndofs) + ierr++; + } } madof[0] = 1; for (size_t n = 1; n < madof.size(); n++) madof[n] += madof[n-1]; - return true; + if (ierr == 0) + return true; + + std::cerr <<" *** SIMbase::addMADOF: Detected "<< ierr <<" nodes with" + <<" conflicting number of DOFs in adjacent patches."<< std::endl; + return false; } @@ -1992,3 +1998,11 @@ const IntVec& SIMbase::getMADOF (unsigned char basis, static IntVec empty; return empty; } + + +void SIMbase::registerDependency (const std::string& name, SIMdependency* sim, + short int nvc, unsigned char basis) +{ + this->registerDependency(sim,name,nvc,myModel, + this->getMADOF(basis,nvc).data()); +} diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index d391b7c11..260492c23 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -561,6 +561,15 @@ class SIMbase : public SIMadmin, public SIMdependency const Vectors& sol, size_t pindx) const; public: + using SIMdependency::registerDependency; + //! \brief Registers a dependency on a field from another SIM object. + //! \param[in] name Name of field we depend on + //! \param[in] sim The SIM object holding the field we depend on + //! \param[in] nvc Number of components in field + //! \param[in] basis The basis the dependent field should be defined on + void registerDependency(const std::string& name, SIMdependency* sim, + short int nvc = 1, unsigned char basis = 1); + //! \brief Extracts all local solution vector(s) for a specified patch. //! \param[in] sol Global primary solution vectors in DOF-order //! \param[in] pindx Local patch index to extract solution vectors for diff --git a/src/SIM/SIMdependency.C b/src/SIM/SIMdependency.C index 6205f9a76..4d75acc98 100644 --- a/src/SIM/SIMdependency.C +++ b/src/SIM/SIMdependency.C @@ -23,13 +23,24 @@ void SIMdependency::registerDependency (SIMdependency* sim, const PatchVec& patches, char diffBasis, int component) { - this->SIMdependency::registerDependency(sim,name,nvc); + this->registerDependency(sim,name,nvc); depFields.back().patches = patches; depFields.back().differentBasis = diffBasis; depFields.back().comp_use = component; } +void SIMdependency::registerDependency (SIMdependency* sim, + const std::string& name, short int nvc, + const PatchVec& patches, + const int* MADOF) +{ + this->registerDependency(sim,name,nvc); + depFields.back().patches = patches; + depFields.back().MADOF = MADOF; +} + + void SIMdependency::registerDependency (SIMdependency* sim, const std::string& name, short int nvc) { @@ -118,46 +129,50 @@ bool SIMdependency::extractPatchDependencies (IntegrandBase* problem, const PatchVec& model, size_t pindx) const { - ASMbase* patch; - DepVector::const_iterator it; - for (it = depFields.begin(); it != depFields.end(); ++it) + for (const Dependency& dp : depFields) { - Vector* lvec = problem->getNamedVector(it->name); + Vector* lvec = problem->getNamedVector(dp.name); if (!lvec) continue; // Ignore fields without corresponding integrand vector - const Vector* gvec = it->sim->getField(it->name); + const Vector* gvec = dp.sim->getField(dp.name); if (!gvec) { std::cerr <<" *** SIMdependency::extractPatchDependencies: \"" - << it->name <<"\" is not a registered field in the simulator \"" - << it->sim->getName() <<"\""<< std::endl; + << dp.name <<"\" is not a registered field in the simulator \"" + << dp.sim->getName() <<"\""<< std::endl; return false; } - else if (gvec->empty()) { + else if (gvec->empty()) + { lvec->clear(); continue; // No error, silently ignore empty fields (treated as zero) } - patch = pindx < it->patches.size() ? it->patches[pindx] : model[pindx]; // See ASMbase::extractNodeVec for interpretation of negative value on basis - int basis = it->components < 0 ? it->components : it->differentBasis; - if (it->differentBasis && it->components != patch->getNoFields(basis)) - patch->extractNodeVec(*gvec,*lvec); + int basis = dp.components < 0 ? dp.components : dp.differentBasis; + ASMbase* pch = pindx < dp.patches.size() ? dp.patches[pindx] : model[pindx]; + if (dp.MADOF) + pch->extractNodalVec(*gvec,*lvec,dp.MADOF); + else if (dp.differentBasis && dp.components != pch->getNoFields(basis)) + pch->extractNodeVec(*gvec,*lvec); else - patch->extractNodeVec(*gvec,*lvec,abs(it->components),basis); - if (it->differentBasis > 0) { - if (it->components == 1) - problem->setNamedField(it->name,Field::create(patch,*lvec, - it->differentBasis, - it->comp_use)); - else - problem->setNamedFields(it->name,Fields::create(patch,*lvec, - it->differentBasis)); - } + pch->extractNodeVec(*gvec,*lvec,abs(dp.components),basis); + #if SP_DEBUG > 2 - std::cout <<"SIMdependency: Dependent field \""<< it->name + std::cout <<"SIMdependency: Dependent field \""<< dp.name <<"\" for patch "<< pindx+1 << *lvec; #endif + if (dp.differentBasis > 0) + { + // Create a field object to handle different interpolation basis + if (dp.components == 1) + problem->setNamedField(dp.name,Field::create(pch,*lvec, + dp.differentBasis, + dp.comp_use)); + else + problem->setNamedFields(dp.name,Fields::create(pch,*lvec, + dp.differentBasis)); + } } return true; diff --git a/src/SIM/SIMdependency.h b/src/SIM/SIMdependency.h index c91e30dd7..47053ae5c 100644 --- a/src/SIM/SIMdependency.h +++ b/src/SIM/SIMdependency.h @@ -42,10 +42,11 @@ class SIMdependency short int components; //!< Number of field components per node short int comp_use; //!< Component to use from field char differentBasis; //!< Toggle usage of an independent basis + const int* MADOF; //!< The MADOF array to associate with field //! \brief Default constructor. Dependency(SIMdependency* s = nullptr, const std::string& f = "", - short int n = 1) : sim(s), name(f), components(n), - comp_use(1), differentBasis(0) {} + short int n = 1) : sim(s), name(f), components(n), comp_use(1), + differentBasis(0), MADOF(nullptr) {} }; //! \brief SIM dependency container @@ -80,6 +81,15 @@ class SIMdependency //! \param[in] sim The SIM object holding the field we depend on //! \param[in] name Name of field we depend on //! \param[in] nvc Number of components in field + //! \param[in] patches The geometry the field is defined over + //! \param[in] MADOF The MADOF array to associate with this field + virtual void registerDependency(SIMdependency* sim, const std::string& name, + short int nvc, const PatchVec& patches, + const int* MADOF); + //! \brief Registers a dependency on a field from another SIM object. + //! \param[in] sim The SIM object holding the field we depend on + //! \param[in] name Name of field we depend on + //! \param[in] nvc Number of components in field virtual void registerDependency(SIMdependency* sim, const std::string& name, short int nvc = 1); diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C index 79ab7c649..75c1d60fc 100644 --- a/src/SIM/SIMoutput.C +++ b/src/SIM/SIMoutput.C @@ -426,15 +426,15 @@ bool SIMoutput::writeGlvV (const Vector& vec, const char* fieldName, IntVec vID; int geomID = myGeomID; - for (size_t i = 0; i < myModel.size(); i++) + for (const ASMbase* pch : myModel) { - if (myModel[i]->empty()) continue; // skip empty patches + if (pch->empty()) continue; // skip empty patches if (msgLevel > 1) - IFEM::cout <<"Writing vector field for patch "<< i+1 << std::endl; + IFEM::cout <<"Writing vector field for patch "<< pch->idx+1 << std::endl; - myModel[i]->extractNodeVec(vec,lovec,ncmp); - if (!myModel[i]->evalSolution(field,lovec,opt.nViz)) + pch->extractNodeVec(vec,lovec,ncmp); + if (!pch->evalSolution(field,lovec,opt.nViz)) return false; if (!myVtf->writeVres(field,++nBlock,++geomID,this->getNoSpaceDim())) @@ -447,6 +447,11 @@ bool SIMoutput::writeGlvV (const Vector& vec, const char* fieldName, } +/*! + This method assumes the scalar field is attached to the first basis + if we are using a mixed basis. +*/ + bool SIMoutput::writeGlvS (const Vector& scl, const char* fieldName, int iStep, int& nBlock, int idBlock) const { @@ -459,16 +464,18 @@ bool SIMoutput::writeGlvS (const Vector& scl, const char* fieldName, Vector lovec; IntVec sID; + int basis = 1; int geomID = myGeomID; - for (size_t i = 0; i < myModel.size(); i++) + for (const ASMbase* pch : myModel) { - if (myModel[i]->empty()) continue; // skip empty patches + if (pch->empty()) continue; // skip empty patches if (msgLevel > 1) - IFEM::cout <<"Writing scalar field for patch "<< i+1 << std::endl; + IFEM::cout <<"Writing scalar field for patch "<< pch->idx+1 << std::endl; - myModel[i]->extractNodeVec(scl,lovec,1); - if (!myModel[i]->evalSolution(field,lovec,opt.nViz)) + int ncmp = scl.size() / this->getNoNodes(basis); + this->extractPatchSolution(scl,lovec,pch,ncmp,basis); + if (!pch->evalSolution(field,lovec,opt.nViz,ncmp)) return false; if (!myVtf->writeNres(field,++nBlock,++geomID)) diff --git a/src/Utility/Test/TestLegendre.C b/src/Utility/Test/TestLegendre.C index 10827d2ba..370c834e5 100644 --- a/src/Utility/Test/TestLegendre.C +++ b/src/Utility/Test/TestLegendre.C @@ -16,6 +16,7 @@ #include "gtest/gtest.h" +#ifdef HAS_BLAS TEST(TestLegendre, GL) { std::vector w, p; @@ -60,6 +61,7 @@ TEST(TestLegendre, GLL) ASSERT_FLOAT_EQ(p[2], 0.447213595499958); ASSERT_FLOAT_EQ(p[3], 1.0); } +#endif TEST(TestLegendre, Eval) @@ -84,14 +86,15 @@ TEST(TestLegendre, Derivative) } +#ifdef HAS_BLAS TEST(TestLegendre, BasisDerivatives) { utl::matrix result1, result2; Legendre::basisDerivatives(3, result1); - const double ref3[3][3] = {{-1.5, 2.0, -0.5}, - {-0.5, 0.0, 0.5}, - { 0.5, -2.0, 1.5}}; + const Real ref3[3][3] = {{-1.5, 2.0, -0.5}, + {-0.5, 0.0, 0.5}, + { 0.5, -2.0, 1.5}}; for (size_t i = 0; i < 3; ++i) for (size_t j = 0; j < 3; ++j) @@ -99,16 +102,16 @@ TEST(TestLegendre, BasisDerivatives) Legendre::basisDerivatives(10, result2); - const double ref10[10][10] = - {{-22.5, 30.43814502928187, -12.17794670742983, 6.943788485133953, - -4.599354761103132, 3.294643033749184, -2.452884175442686, 1.829563931903247, + const Real ref10[10][10] = + {{-22.5, 30.43814502928187, -12.17794670742983, 6.943788485133953, + -4.599354761103132, 3.294643033749184, -2.452884175442686, 1.829563931903247, -1.275954836092663, 0.5}, {-5.074064702978062, 0.0, 7.185502869705838, -3.351663862746774, 2.078207994036418, -1.444948448751456, 1.059154463645442, -0.783239293137909, 0.5437537382357057, -0.2127027580091889}, { 1.203351992852207, -4.259297354965217, 0.0, 4.368674557010181, -2.104350179413155, 1.334915483878251, -0.9366032131394465, - 0.6767970871960859, -0.4642749589081571, 0.1807865854892499}, + 0.6767970871960859, -0.4642749589081571, 0.1807865854892499}, {-0.528369376820273, 1.529902638181603, -3.364125868297819, 0.0, 3.387318101202445, -1.64649408398706, 1.046189365502494, -0.7212373127216044, 0.4834623263339481, -0.1866457893937361}, @@ -125,7 +128,7 @@ TEST(TestLegendre, BasisDerivatives) 0.9366032131394465, -1.334915483878251, 2.104350179413155, -4.368674557010181, 0.0, 4.259297354965217, -1.203351992852207}, - {0.2127027580091889, -0.5437537382357057, 0.783239293137909, + { 0.2127027580091889, -0.5437537382357057, 0.783239293137909, -1.059154463645442, 1.444948448751456, -2.078207994036418, 3.351663862746774, -7.185502869705838, 0.0, 5.074064702978062}, {-0.5, 1.275954836092663, -1.829563931903247, @@ -136,3 +139,4 @@ TEST(TestLegendre, BasisDerivatives) for (size_t j = 0; j < 10; ++j) EXPECT_FLOAT_EQ(result2(i+1, j+1), ref10[i][j]); } +#endif