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