1 /** @file gsVisitorMass.h 2 3 @brief Mass visitor for assembling element mass matrix 4 5 This file is part of the G+Smo library. 6 7 This Source Code Form is subject to the terms of the Mozilla Public 8 License, v. 2.0. If a copy of the MPL was not distributed with this 9 file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 Author(s): A. Mantzaflaris 12 */ 13 14 #pragma once 15 16 namespace gismo 17 { 18 19 /** @brief The visitor computes element mass integrals 20 * 21 * Assembles the bilinear term 22 * \f[ (u,v)_\Omega, \f] 23 * where \f$u\f$ is the trial function and \f$v\f$ is the test function. 24 * 25 * @ingroup Assembler 26 */ 27 28 template <class T> 29 class gsVisitorMass 30 { 31 public: 32 33 /// Constructor gsVisitorMass()34 gsVisitorMass() 35 { } 36 37 /// @brief Constructor 38 /// 39 /// @param pde Reference to \a gsPde object (is ignored) gsVisitorMass(const gsPde<T> & pde)40 gsVisitorMass(const gsPde<T> & pde) 41 { GISMO_UNUSED(pde); } 42 43 /// Initialize initialize(const gsBasis<T> & basis,const index_t,const gsOptionList & options,gsQuadRule<T> & rule)44 void initialize(const gsBasis<T> & basis, 45 const index_t , 46 const gsOptionList & options, 47 gsQuadRule<T> & rule) 48 { 49 // Setup Quadrature (harmless slicing occurs) 50 rule = gsQuadrature::get(basis, options); // harmless slicing occurs here 51 52 // Set Geometry evaluation flags 53 md.flags = NEED_MEASURE; 54 } 55 56 /// Evaluate on element evaluate(const gsBasis<T> & basis,const gsGeometry<T> & geo,gsMatrix<T> & quNodes)57 inline void evaluate(const gsBasis<T> & basis, // to do: more unknowns 58 const gsGeometry<T> & geo, 59 // todo: add element here for efficiency 60 gsMatrix<T> & quNodes) 61 { 62 md.points = quNodes; 63 // Compute the active basis functions 64 // Assumes actives are the same for all quadrature points on the current element 65 basis.active_into(md.points.col(0), actives); 66 const index_t numActive = actives.rows(); 67 68 // Evaluate basis functions on element 69 basis.eval_into(md.points, basisData); 70 71 // Compute geometry related values 72 geo.computeMap(md); 73 74 // Initialize local matrix/rhs 75 localMat.setZero(numActive, numActive); 76 } 77 78 /// Assemble on element assemble(gsDomainIterator<T> &,gsVector<T> const & quWeights)79 inline void assemble(gsDomainIterator<T> & , 80 gsVector<T> const & quWeights) 81 { 82 localMat.noalias() = 83 basisData * quWeights.asDiagonal() * 84 md.measures.asDiagonal() * basisData.transpose(); 85 } 86 87 /// Adds the contributions to the sparse system localToGlobal(const index_t patchIndex,const std::vector<gsMatrix<T>> & eliminatedDofs,gsSparseSystem<T> & system)88 inline void localToGlobal(const index_t patchIndex, 89 const std::vector<gsMatrix<T> > & eliminatedDofs, 90 gsSparseSystem<T> & system) 91 { 92 // Map patch-local DoFs to global DoFs 93 system.mapColIndices(actives, patchIndex, actives); 94 95 // Add contributions to the system matrix 96 system.pushToMatrix(localMat, actives, eliminatedDofs.front(), 0, 0); 97 } 98 99 /* ----------------------- to be removed later*/ 100 initialize(const gsBasis<T> & basis,gsQuadRule<T> & rule)101 void initialize(const gsBasis<T> & basis, 102 gsQuadRule<T> & rule) 103 { 104 gsVector<short_t> numQuadNodes( basis.dim() ); 105 for (short_t i = 0; i < basis.dim(); ++i) 106 numQuadNodes[i] = basis.degree(i) + 1; 107 108 // Setup Quadrature 109 rule = gsGaussRule<T>(numQuadNodes);// harmless slicing occurs here 110 111 // Set Geometry evaluation flags 112 md.flags = NEED_MEASURE; 113 } 114 115 localToGlobal(const gsDofMapper & mapper,const gsMatrix<T> & eliminatedDofs,const index_t patchIndex,gsSparseMatrix<T> & sysMatrix,gsMatrix<T> & rhsMatrix)116 void localToGlobal(const gsDofMapper & mapper, 117 const gsMatrix<T> & eliminatedDofs, 118 const index_t patchIndex, 119 gsSparseMatrix<T> & sysMatrix, 120 gsMatrix<T> & rhsMatrix ) 121 { 122 mapper.localToGlobal(actives, patchIndex, actives); 123 124 const index_t numActive = actives.rows(); 125 126 for (index_t i = 0; i < numActive; ++i) 127 { 128 const int ii = actives(i,0); // N_i 129 130 if ( mapper.is_free_index(ii) ) 131 { 132 for (index_t j = 0; j < numActive; ++j) 133 { 134 const int jj = actives(j,0); // N_j 135 if ( mapper.is_free_index(jj) ) 136 //if ( jj <= ii ) // only store lower triangular part 137 sysMatrix.coeffRef(ii, jj) += localMat(i, j); // N_i*N_j 138 } 139 } 140 } 141 } 142 143 144 protected: 145 146 // Basis values 147 gsMatrix<T> basisData; 148 gsMatrix<index_t> actives; 149 150 // Local matrix 151 gsMatrix<T> localMat; 152 153 gsMapData<T> md; 154 }; 155 156 157 } // namespace gismo 158