| Directory: | ./ |
|---|---|
| File: | FiniteElement/elementMatrix.hpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 22 | 22 | 100.0% |
| Branches: | 8 | 16 | 50.0% |
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // ______ ______ _ _ _____ ______ | ||
| 2 | // | ____| ____| | (_)/ ____| | ____| | ||
| 3 | // | |__ | |__ | | _| (___ ___| |__ | ||
| 4 | // | __| | __| | | | |\___ \ / __| __| | ||
| 5 | // | | | |____| |____| |____) | (__| |____ | ||
| 6 | // |_| |______|______|_|_____/ \___|______| | ||
| 7 | // Finite Elements for Life Sciences and Engineering | ||
| 8 | // | ||
| 9 | // License: LGL2.1 License | ||
| 10 | // FELiScE default license: LICENSE in root folder | ||
| 11 | // | ||
| 12 | // Main authors: J. Foulon & J-F. Gerbeau & V. Martin & Miguel Fernandez | ||
| 13 | // | ||
| 14 | |||
| 15 | #ifndef _ELEMMAT_H_INCLUDED | ||
| 16 | #define _ELEMMAT_H_INCLUDED | ||
| 17 | |||
| 18 | // System includes | ||
| 19 | #include <vector> | ||
| 20 | #include <ostream> | ||
| 21 | #include <atomic> | ||
| 22 | |||
| 23 | // External includes | ||
| 24 | |||
| 25 | // Project includes | ||
| 26 | #include "Core/felisce.hpp" | ||
| 27 | #include "Core/shared_pointers.hpp" | ||
| 28 | #include "FiniteElement/elementField.hpp" | ||
| 29 | #include "FiniteElement/elementVector.hpp" | ||
| 30 | #include "FiniteElement/curBaseFiniteElement.hpp" | ||
| 31 | |||
| 32 | namespace felisce | ||
| 33 | { | ||
| 34 | ///@name felisce Globals | ||
| 35 | ///@{ | ||
| 36 | |||
| 37 | ///@} | ||
| 38 | ///@name Type Definitions | ||
| 39 | ///@{ | ||
| 40 | |||
| 41 | ///@} | ||
| 42 | ///@name Enum's | ||
| 43 | ///@{ | ||
| 44 | |||
| 45 | ///@} | ||
| 46 | ///@name Functions | ||
| 47 | ///@{ | ||
| 48 | |||
| 49 | ///@} | ||
| 50 | ///@name felisce Classes | ||
| 51 | ///@{ | ||
| 52 | class ElementMatrix | ||
| 53 | { | ||
| 54 | public: | ||
| 55 | ///@name Type Definitions | ||
| 56 | ///@{ | ||
| 57 | |||
| 58 | /// Pointer definition of ElementMatrix | ||
| 59 | FELISCE_CLASS_POINTER_DEFINITION(ElementMatrix); | ||
| 60 | |||
| 61 | ///@} | ||
| 62 | ///@name Life Cycle | ||
| 63 | ///@{ | ||
| 64 | |||
| 65 | // default constructor | ||
| 66 |
6/12✓ Branch 6 taken 120 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 120 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 120 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 120 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 120 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 120 times.
✗ Branch 22 not taken.
|
120 | ElementMatrix() = default; |
| 67 | |||
| 68 | //! constructor for an arbitrary number of finite elements | ||
| 69 | ElementMatrix(const std::vector<const CurBaseFiniteElement*>& fe, const std::vector<std::size_t>& nbr, const std::vector<std::size_t>& nbc); | ||
| 70 | |||
| 71 | //! constructor for an arbitrary number of finite elements | ||
| 72 | ElementMatrix(const std::vector<const CurBaseFiniteElement*>& fer, const std::vector<std::size_t>& nbr, const std::vector<const CurBaseFiniteElement*>& fec, const std::vector<std::size_t>& nbc); | ||
| 73 | |||
| 74 | // constructor for 1 finite element | ||
| 75 | FELISCE_DEPRECATED_MESSAGE("WARNING:: This constructor will be removed in the near future. Please use the main constructor which consider constant vector.") | ||
| 76 | ElementMatrix(const CurBaseFiniteElement& fe, const std::size_t nbr1, const std::size_t nbc1); | ||
| 77 | |||
| 78 | // constructor for 2 finite elements | ||
| 79 | FELISCE_DEPRECATED_MESSAGE("WARNING:: This constructor will be removed in the near future. Please use the main constructor which consider constant vector.") | ||
| 80 | ElementMatrix(const CurBaseFiniteElement& fe1, const std::size_t nbr1, const std::size_t nbc1, | ||
| 81 | const CurBaseFiniteElement& fe2, const std::size_t nbr2, const std::size_t nbc2); | ||
| 82 | |||
| 83 | // constructor for 3 finite elements | ||
| 84 | FELISCE_DEPRECATED_MESSAGE("WARNING:: This constructor will be removed in the near future. Please use the main constructor which consider constant vector.") | ||
| 85 | ElementMatrix(const CurBaseFiniteElement& fe1, const std::size_t nbr1, const std::size_t nbc1, | ||
| 86 | const CurBaseFiniteElement& fe2, const std::size_t nbr2, const std::size_t nbc2, | ||
| 87 | const CurBaseFiniteElement& fe3, const std::size_t nbr3, const std::size_t nbc3 ); | ||
| 88 | |||
| 89 | |||
| 90 | //! Copy constructor, to avoid warning due to user-declared destructor. | ||
| 91 | ElementMatrix(const ElementMatrix& rOther); | ||
| 92 | |||
| 93 | // default destructor | ||
| 94 | 72183 | ~ElementMatrix() = default; | |
| 95 | |||
| 96 | ///@} | ||
| 97 | ///@name Operators | ||
| 98 | ///@{ | ||
| 99 | |||
| 100 | //! Assignment, to avoid warning due to user-declared destructor. | ||
| 101 | ElementMatrix& operator=(const ElementMatrix& rOther); | ||
| 102 | |||
| 103 | /*! | ||
| 104 | * \brief operator += | ||
| 105 | * | ||
| 106 | * This operator acts upon the underlying UBlasMatrix and is is assumed all the other attributes are the same (it | ||
| 107 | * is checked by several assert in debug mode). | ||
| 108 | */ | ||
| 109 | void operator+=(const ElementMatrix& mat); | ||
| 110 | |||
| 111 | /*! | ||
| 112 | * \brief operator -= | ||
| 113 | * | ||
| 114 | * This operator acts upon the underlying UBlasMatrix and is is assumed all the other attributes are the same (it | ||
| 115 | * is checked by several assert in debug mode). | ||
| 116 | */ | ||
| 117 | void operator-=(const ElementMatrix& mat); | ||
| 118 | |||
| 119 | //! Operator multiplication by a scalar. | ||
| 120 | void operator*=(const double factor); | ||
| 121 | |||
| 122 | /** | ||
| 123 | * @brief Access operators | ||
| 124 | */ | ||
| 125 | const double& operator () (const std::size_t i, const std::size_t j) const { | ||
| 126 | return m_mat(i, j); | ||
| 127 | } | ||
| 128 | |||
| 129 | 2079 | double& operator () (const std::size_t i, const std::size_t j) { | |
| 130 | 2079 | return m_mat(i, j); | |
| 131 | } | ||
| 132 | |||
| 133 | ///@} | ||
| 134 | ///@name Operations | ||
| 135 | ///@{ | ||
| 136 | |||
| 137 | void duplicateFrom(const ElementMatrix& mat); | ||
| 138 | |||
| 139 | 341183149 | UBlasMatrixRange matBlock( std::size_t i, std::size_t j ) { | |
| 140 | 341183149 | return getMatBlock(m_mat, i, j); | |
| 141 | } | ||
| 142 | |||
| 143 | 437464137 | UBlasMatrixRange getMatBlock(UBlasMatrix& mat, std::size_t i, std::size_t j ) { | |
| 144 |
1/2✓ Branch 4 taken 437464137 times.
✗ Branch 5 not taken.
|
437464137 | return UBlasMatrixRange(mat,UBlasRange( m_firstRow[i], m_firstRow[i]+m_numRow[i]), |
| 145 |
1/2✓ Branch 5 taken 437464137 times.
✗ Branch 6 not taken.
|
1312392411 | UBlasRange( m_firstCol[ j ],m_firstCol[ j ]+m_numCol[ j ])); |
| 146 | } | ||
| 147 | |||
| 148 | /** | ||
| 149 | * @brief Size 1 | ||
| 150 | */ | ||
| 151 | 238 | std::size_t size1() const { | |
| 152 | 238 | return m_mat.size1(); | |
| 153 | } | ||
| 154 | |||
| 155 | /** | ||
| 156 | * @brief Size 2 | ||
| 157 | */ | ||
| 158 | 2030 | std::size_t size2() const { | |
| 159 | 2030 | return m_mat.size2(); | |
| 160 | } | ||
| 161 | |||
| 162 | /** | ||
| 163 | * @brief Sets to zero the main matrix | ||
| 164 | */ | ||
| 165 | 54995499 | void zero() { | |
| 166 | 54995499 | m_mat.clear(); | |
| 167 | 54995499 | } | |
| 168 | |||
| 169 | //Operator div(grad symmetric): div(\eps(\u)) with \eps(\u) = 1/2 (grad u + grad^T u) | ||
| 170 | // \int \eps(\phi_i):\eps(\phi_j) | ||
| 171 | void eps_phi_i_eps_phi_j(const double coef, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 172 | |||
| 173 | //Operator div(grad symmetric): div(a \eps(\u)) with \eps(\u) = 1/2 (grad u + grad^T u) | ||
| 174 | // \int a \eps(\phi_i):\eps(\phi_j) | ||
| 175 | void a_eps_phi_i_eps_phi_j(const double coef,const ElementField& scalfct,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 176 | |||
| 177 | void grad_phi_i_grad_phi_j(const double coef,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 178 | |||
| 179 | void grad_phi_i_grad_phi_j(const double coef,const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 180 | |||
| 181 | void phi_j_grad_psi_i(const double coef, const CurrentFiniteElement& fei,const CurrentFiniteElement& fej,const std::size_t iblock, const std::size_t jblock); | ||
| 182 | |||
| 183 | void phi_i_grad_psi_j(const double coef, const CurrentFiniteElement& fei,const CurrentFiniteElement& fej,const std::size_t iblock, const std::size_t jblock); | ||
| 184 | |||
| 185 | void a_grad_phi_i_grad_phi_j(const double coef,const ElementField& scalfct,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 186 | |||
| 187 | void a_grad_phi_i_grad_phi_j(const double coef,const ElementField& scalfct,const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 188 | |||
| 189 | //--------------------------------------------------------------------- | ||
| 190 | // compute \int_K \grad \phi_i \cdot \vec \phi_j | ||
| 191 | void grad_phi_i_dot_vec_phi_j(const double coef, const ElementField& vec, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 192 | |||
| 193 | //--------------------------------------------------------------------- | ||
| 194 | // compute \int_K eps \phi_i \cdot \vec \phi_j | ||
| 195 | void eps_phi_i_dot_vec_phi_j(const double coef, const ElementField& vec, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 196 | |||
| 197 | //--------------------------------------------------------------------- | ||
| 198 | // compute \int_K \grad \phi_j \cdot \vec \phi_i | ||
| 199 | void phi_i_grad_phi_j_dot_vec(const double coef, const ElementField& vec, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 200 | |||
| 201 | //--------------------------------------------------------------------- | ||
| 202 | // compute \int_K \eps \phi_j \cdot \vec \phi_i | ||
| 203 | void phi_i_eps_phi_j_dot_vec(const double coef, const ElementField& vec, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 204 | |||
| 205 | //--------------------------------------------------------------------- | ||
| 206 | // compute \int_K \eps \phi_j \vec \cdot \eps \phi_i \vec | ||
| 207 | void eps_phi_j_vec_dot_eps_phi_i_vec(const double coef, const ElementField& vec, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 208 | |||
| 209 | //--------------------------------------------------------------------- | ||
| 210 | // compute \int_K \grad \phi_j \vec \cdot \grad \phi_i \vec | ||
| 211 | void grad_phi_j_vec_dot_grad_phi_i_vec(const double coef, const ElementField& vec, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 212 | |||
| 213 | //---------------------------------------------------------------------- | ||
| 214 | // TODO: instead of diffTensor we should pass some kind of element field | ||
| 215 | // diffTensor is in practice a vector dof_field. | ||
| 216 | void tau_grad_phi_i_tau_grad_phi_j(const double coef,const std::vector<double>& diffTensor,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num,std::size_t dim = 3); | ||
| 217 | |||
| 218 | void tau_grad_phi_i_tau_grad_phi_j(const double coef,const std::vector<double>& diffTensor,const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 219 | |||
| 220 | void tau_orthotau_grad_phi_i_grad_phi_j(const double coef,const std::vector<double>& diffTensor,const std::vector<double>& angle,const CurvilinearFiniteElement& fe,const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 221 | |||
| 222 | //! Operator used in linear elasticity problem | ||
| 223 | void grad_phi_i_GradientBasedElastTensor_grad_phi_j(const double coef,const CurrentFiniteElement& fe, UBlasMatrix& GradientBasedElastTensor, const std::size_t iblock, const std::size_t jblock); | ||
| 224 | |||
| 225 | /*! | ||
| 226 | * \brief Operator used in hyperelasticity problem. | ||
| 227 | */ | ||
| 228 | void grad_phi_i_GradientBasedHyperElastTensor_grad_phi_j(const double coef, | ||
| 229 | const CurrentFiniteElement& feDisplacement, | ||
| 230 | int quadraturePointIndex, | ||
| 231 | UBlasMatrix& GradientBasedHyperElastTensor, const std::size_t iblock, const std::size_t jblock); | ||
| 232 | |||
| 233 | void phi_i_phi_j(const double coef,const CurBaseFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 234 | |||
| 235 | void setZerosLung(const double coef,const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 236 | |||
| 237 | void phi_i_phi_j_lumped(const double coef,const CurBaseFiniteElement& fe,const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 238 | |||
| 239 | void a_phi_i_phi_j(const double coef,const ElementField& scalfct,const CurBaseFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 240 | |||
| 241 | void penalty_contact_phi_i_scalar_n_phi_j_scalar_n(const double coef,const std::vector<double>& normal, const CurvilinearFiniteElement& fe,const ElementField& disp, const ElementField& gap); | ||
| 242 | |||
| 243 | void penalty_phi_i_scalar_n_phi_j_scalar_n(const double coef, CurvilinearFiniteElement& fe, const ElementField& vel, std::size_t iblock, std::size_t jblock); | ||
| 244 | |||
| 245 | void nitsche_contact_phi_i_scalar_n_phi_j_scalar_n(const double coef,const std::vector<double>& normal, const CurvilinearFiniteElement& fe,const ElementField& disp, const ElementField& gap, const ElementField& lambda); | ||
| 246 | |||
| 247 | void nitsche_contact_phi_i_scalar_n_phi_j_scalar_n(const double coef, const CurvilinearFiniteElement& fe,const ElementField& disp, const ElementField& gap, const ElementField& normal, const ElementField& lambda); | ||
| 248 | |||
| 249 | void phi_i_psi_j_scalar_n(const double coef,const CurvilinearFiniteElement& fe1, const CurvilinearFiniteElement& fe2, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 250 | |||
| 251 | void u_grad_phi_j_phi_i(const double coef,const ElementField& vel,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num, bool transpose=false); | ||
| 252 | |||
| 253 | void u_grad_phi_j_phi_i(const double coef,const ElementField& vel, const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 254 | |||
| 255 | void u_grad_phi_i_grad_psi_j(const double coef,const ElementField& vel,const CurrentFiniteElement& fe1, const CurrentFiniteElement& fe2, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 256 | |||
| 257 | //---------------------------------------------------------------------- | ||
| 258 | // Operator \intOmega coef (rot(u^n+1) x u^n) \cdot v - coef/2 (u^n+1 \cdot u^n) div(v) | ||
| 259 | // Used in linearProblemNS : \int_{\Omega} \rho rot(u^{n+1} x u^n \cdot v - \rho/2 u^{n+1} \cdot u^n) div(v) | ||
| 260 | // with coef = \rho | ||
| 261 | void convective_rot(const double coef, const ElementField& vel, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 262 | |||
| 263 | //---------------------------------------------------------------------- | ||
| 264 | // operator \int_{Sigma} u^{n+1} Grad u^{n} \dot v | ||
| 265 | void phi_grad_u(const double coef,const ElementField& vel,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 266 | |||
| 267 | void phi_grad_u_minus_grad_u_phi(const double coef,const ElementField& vel,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 268 | |||
| 269 | // Temam's operator | ||
| 270 | // Operator \intOmega (\nabla \cdot u^n) u^n+1 v (Temam's operator to stabilize the convective term) | ||
| 271 | void div_u_phi_j_phi_i(const double coef, const ElementField& vel, const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 272 | |||
| 273 | //! Operators \f$ \int_{\Sigma} (\grad \phi_i \dot n) \psi_j \f$, (and its transpose) | ||
| 274 | //! where $\Sigma$ is a boundary element. (For Nitsche's treatment...) | ||
| 275 | void grad_phi_i_dot_n_psi_j(const CurrentFiniteElementWithBd& fewbd1, const CurrentFiniteElementWithBd& fewbd2, int iblockBd, int numblockBd, const std::size_t iblock, const std::size_t jblock,const std::size_t num, const double coef1,const double coef2); | ||
| 276 | |||
| 277 | //! Operator \f$ \int_{\Sigma} (\grad \phi_i \dot n) (\grad \phi_j \dot n) \f$, | ||
| 278 | //! where $\Sigma$ is a boundary element. (For Nitsche's treatment...) | ||
| 279 | void grad_phi_i_dot_n_grad_phi_j_dot_n(const double coef,const CurrentFiniteElementWithBd& fe, int iblockBd, int numblockBd, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 280 | |||
| 281 | //! Operator \int_{\Sigma} (\grad \phi_i \dot n) (\grad \phi_j \dot n) | ||
| 282 | //! Use to compute the ghost penalty stabilization term | ||
| 283 | //! where $\Sigma$ is a boundary element. | ||
| 284 | //! iBdPhi is the local id of the boundary for the element fewbdphi | ||
| 285 | //! iBdPsi is the local id of the same boundary but for the element fewbdpsi | ||
| 286 | void grad_phi_j_dot_n_grad_psi_i_dot_n(const double coef, const CurrentFiniteElementWithBd& fewbdphi, const CurrentFiniteElementWithBd& fewbdpsi, felInt iBdPhi, felInt iBdPsi, const std::size_t iblock, const std::size_t jblock, const std::size_t numBlock); | ||
| 287 | |||
| 288 | //! Operator \int_{\Sigma} (\grad \phi_i) (\grad \phi_j) | ||
| 289 | //! Use to compute the pressure jump for face-oriented stabilization | ||
| 290 | //! where $\Sigma$ is a boundary element. | ||
| 291 | //! iBdPhi is the local id of the boundary for the element fewbdphi | ||
| 292 | //! iBdPsi is the local id of the same boundary but for the element fewbdpsi | ||
| 293 | void grad_phi_j_grad_psi_i(const double coef, const CurrentFiniteElementWithBd& fewbdphi, const CurrentFiniteElementWithBd& fewbdpsi, felInt iBdPhi, felInt iBdPsi, const std::size_t iblock, const std::size_t jblock, const std::size_t numBlock); | ||
| 294 | |||
| 295 | void div_phi_j_div_phi_i(const double coef,const CurrentFiniteElement& fe, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 296 | |||
| 297 | /*! | ||
| 298 | Operator \f$ \int \psi_j \mbox{div} \phi_i \f$, multiplied by coef1 in blocks (iblock,jblock), ..., (iblock+numCoor,jblock), | ||
| 299 | and \f$ \int \psi_i \mbox{div} \phi_j \f$ calculated by transposing the first term, multiplied by coef2 in blocks (jblock,iblock),...,(jblock,iblock+numCoor) | ||
| 300 | Example : \int p \nabla \cdot v and \int q \nabla \cdot u in Stokes problem | ||
| 301 | */ | ||
| 302 | void psi_j_div_phi_i_and_psi_i_div_phi_j (const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2, const std::size_t iblock, const std::size_t jblock,const double coef1,const double coef2); | ||
| 303 | /* | ||
| 304 | Operator \f$ coeff \int \psi_j \nabla \cdot \phi_i \f$ | ||
| 305 | */ | ||
| 306 | void psi_j_div_phi_i(const double coef, const CurrentFiniteElement& fei,const CurrentFiniteElement& fej,const std::size_t iblock, const std::size_t jblock); | ||
| 307 | /* | ||
| 308 | Operator \f$ coeff \int \psi_i \nabla \cdot \phi_j \f$ | ||
| 309 | */ | ||
| 310 | void psi_i_div_phi_j(const double coef, const CurrentFiniteElement& fej,const CurrentFiniteElement& fei,const std::size_t iblock, const std::size_t jblock); | ||
| 311 | |||
| 312 | //--------------------------------------------------------------------------------------------------------------- | ||
| 313 | // Some explanation about SUPG: | ||
| 314 | // noting T:time, L:length, M:mass | ||
| 315 | // WARNING: tau has a dimension of T*L^3*M^{-1} -> remember this when you use P1-P1 | ||
| 316 | // It's advised to choose your units knowing this. | ||
| 317 | //--------------------------------------------------------------------------------------------------------------- | ||
| 318 | void stab_supg(const double dt, const double stab1, const double stab2, const double density, const double visc,const CurrentFiniteElement& fe, const ElementField& vel,const ElementField& rhs,ElementVector& elemVec, double& Re_elem,double& tau, const int type=1); | ||
| 319 | |||
| 320 | void stab_supgNSHeat(const double dt, const double stab1, const double alpha, const double density, const double visc, const CurrentFiniteElement& fe, const ElementField& vel, const int type,const std::vector<double> & gravity); | ||
| 321 | |||
| 322 | void stab_supgAdvDiffCT(const double dt, const double stab1, const double stab2, const double density, const double visc, const CurrentFiniteElement& fe,const ElementField& vel); | ||
| 323 | |||
| 324 | void stab_supgProjCT(const double dt, const double stab1, const double density, const double visc, const CurrentFiniteElement& fe,const ElementField& vel); | ||
| 325 | |||
| 326 | void stab_supgMatrixFS(const double dt, const double stabSUPG, const double density, const double viscosity, const double c0, const double c1, const double c2, const CurrentFiniteElement& fe, const ElementField& vel); | ||
| 327 | |||
| 328 | // outflow stabilization a la Bertoglio et al. (in progress) | ||
| 329 | void f_dot_n_grad_phi_i_grad_phi_j(const double coef,const CurvilinearFiniteElement& fe, const ElementField& elfield,const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 330 | |||
| 331 | //void outflow_stabilization(const double coef, const CurvilinearFiniteElement& fe, const ElementField& elfield, const std::size_t iblock, const std::size_t jblock,const std::size_t num); | ||
| 332 | // Operator: int_Sigma [ ( u^{n} \cdot normal ) ( u^{n+1}*v ) ] | ||
| 333 | // If you want to apply outflow stabilization method as the following example: | ||
| 334 | // Example : if f_dot_n = u^{n} * normal < 0, | ||
| 335 | // adding - (rho/2) * int_Sigma [ ( u^{n} \cdot normal ) ( u^{n+1}*v ) ] | ||
| 336 | // You should activate stabilization flag in your data file (see Core/felisceParam.cpp -> addedBoundaryFlag | ||
| 337 | void f_dot_n_phi_i_phi_j(const double coef, const CurvilinearFiniteElement& fe, const ElementField& elfield, const std::size_t iblock, const std::size_t jblock,const std::size_t num, int stabFlag = 0); | ||
| 338 | |||
| 339 | // inflow stabilization for Nitsche-XFEM formulation | ||
| 340 | // Same as f_dot_n_phi_i_phi_j but with currentFiniteElement instead of CurvilinearFiniteElement | ||
| 341 | // Used in Nitsche-XFEM formulation | ||
| 342 | // elt_f is the function f | ||
| 343 | // elt_n is the "normal". This is a constant elementField with as many components as elt_f | ||
| 344 | void f_dot_n_phi_i_phi_j(const double coef, const CurrentFiniteElement& fe, const ElementField& elt_f, const ElementField& elt_n, const std::size_t iblock, const std::size_t jblock, const std::size_t num, const felInt stabFlag); | ||
| 345 | |||
| 346 | /*! | ||
| 347 | Operator | ||
| 348 | \int_{\Sigma} ( \|\phi \dot n \| - (\phi \dot n) ) \phi_i \phi_j | ||
| 349 | for selecting the inflow boundary (where u.n < 0) | ||
| 350 | */ | ||
| 351 | void abs_u_dot_n_phi_i_phi_j(const double coef, ElementField& vel, CurvilinearFiniteElement& febd, const std::size_t iblock, const std::size_t jblock, const int num); | ||
| 352 | |||
| 353 | // Operator: coef * int_Sigma [ ( u \times normale ) ( v \times normale ) ] (in 3D) | ||
| 354 | // To impose the tangential velocity to be zero (for now) | ||
| 355 | // The labels must be specified in the datafile (see Core/felisceParam.cpp -> compTangBClabel ) | ||
| 356 | void phi_i_times_n_phi_j_times_n(const double coef, const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const int numComp); | ||
| 357 | |||
| 358 | // Operator: coef * int_Sigma [ ( u \cdot normale ) ( v \cdot normale ) ] | ||
| 359 | void phi_i_dot_n_phi_j_dot_n(const double coef, const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const int numComp); | ||
| 360 | |||
| 361 | // Operator: coef * int_Sigma [ ( u \cdot tangent ) ( v \cdot tangent ) ] | ||
| 362 | void phi_i_dot_t_phi_j_dot_t(const double coef, const CurvilinearFiniteElement& fe, const std::size_t iblock, const std::size_t jblock, const int numComp); | ||
| 363 | |||
| 364 | // int_Sigma [coeff eps(phi_j)n dot n phi_i dot n ] (i=test function) | ||
| 365 | void eps_phi_j_n_dot_n_phi_i_dot_n(const double coeff, const CurrentFiniteElementWithBd& fe,int iblockBd,int numblockBd,const std::size_t iblock,const std::size_t jblock, const std::size_t num); | ||
| 366 | |||
| 367 | // int_Sigma [coeff psi_j phi_i dot n ] (i=test function) something similar already there in phi_i_psi_j_scalar_n | ||
| 368 | void psi_j_phi_i_dot_n(const double coeff,const CurvilinearFiniteElement& feVel, const CurvilinearFiniteElement& fePres,const std::size_t iblock,std::size_t presBlock, const std::size_t num); | ||
| 369 | |||
| 370 | //for each node the std::vector P*N is defined, then it is interpolated via the FE basis into the element | ||
| 371 | void psi_j_phi_i_dot_n_on_nodes(const double coeff,const ElementField& normal, const CurvilinearFiniteElement& feVel, const CurvilinearFiniteElement& fePres,const std::size_t iblock, const std::size_t presBlock); | ||
| 372 | |||
| 373 | void sGrad_psi_j_tensor_sGrad_phi_i(const double coef, const ElementField& normalField, const CurvilinearFiniteElement& fe, const UBlasMatrix& tensor); | ||
| 374 | |||
| 375 | /** | ||
| 376 | * @brief surface laplacian operator for scalar quantities | ||
| 377 | * @param[in] coef coefficient in front of the integral | ||
| 378 | * @param[in] fe finite element of the surface | ||
| 379 | * @param[in] tensor the tensor between the covariant gradients: if you put the inverse of the first fundamental form you obtain a laplacian operator. | ||
| 380 | * @param[in] blockId block where to put the result. | ||
| 381 | */ | ||
| 382 | void sGrad_psi_j_tensor_sGrad_phi_i_for_scalar(const double coef, const CurvilinearFiniteElement& fe, const UBlasMatrix& tensor, const std::size_t blockId); | ||
| 383 | |||
| 384 | void sGrad_psi_j_tensor_sGrad_phi_i(const double coef, const ElementField& elemCoef, const ElementField& normalField, const CurvilinearFiniteElement& fe, const UBlasMatrix& tensor); | ||
| 385 | |||
| 386 | // used in the choroid | ||
| 387 | void dpsi_j_dh_psi_i(const double coef, const CurrentFiniteElement& fei,const CurrentFiniteElement& fej,std::size_t h, const std::size_t iblock, const std::size_t jblock); | ||
| 388 | |||
| 389 | void transpirationSurfStab( const double coef, const ElementField& elemCoef, const ElementField& normalField, const CurvilinearFiniteElement& fe, const UBlasMatrix& tensor); | ||
| 390 | |||
| 391 | // feVel1 correspond to the Volume element, while feVel2 to the surface element | ||
| 392 | // this is used to create a rectangular elementary matrix | ||
| 393 | ElementMatrix(const CurBaseFiniteElement& feVel1, const CurBaseFiniteElement& feVel2, std::size_t nbr1, std::size_t nbc1, | ||
| 394 | const CurBaseFiniteElement& fePres, std::size_t nbr2, std::size_t nbc2); | ||
| 395 | |||
| 396 | void followerPressure(const double coef, const ElementField& press, const ElementField& disp, const CurvilinearFiniteElement& fe); | ||
| 397 | |||
| 398 | void velocityConstraint(const double coef, const ElementField& disp, const CurvilinearFiniteElement& fe, int iblock); | ||
| 399 | |||
| 400 | //---------------------------------------------------------------------- | ||
| 401 | //! Operator for DG on faces \int_{\Sigma} \psi_j \phi_i $, | ||
| 402 | //! where $\Sigma$ is a boundary element. | ||
| 403 | //---------------------------------------------------------------------- | ||
| 404 | void psi_j_phi_i(const double coef, const CurrentFiniteElement& fePhi, const CurrentFiniteElement& fePsi, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 405 | |||
| 406 | //---------------------------------------------------------------------- | ||
| 407 | //! Operator for DG on faces $ \int_{\Sigma} psi_j_eps_phi_i_dot_vec $, | ||
| 408 | //! where $\Sigma$ is a boundary element. | ||
| 409 | //--------------------------------------------------------------------- | ||
| 410 | void psi_j_eps_phi_i_dot_vec(const double coef, const ElementField& vec, const CurrentFiniteElement& fePhi, const CurrentFiniteElement& fePsi, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 411 | |||
| 412 | //---------------------------------------------------------------------- | ||
| 413 | //! Operator for DG on faces $ \int_{\Sigma} psi_j_grad_phi_i_dot_vec $, | ||
| 414 | //! where $\Sigma$ is a boundary element. | ||
| 415 | //--------------------------------------------------------------------- | ||
| 416 | void psi_j_grad_phi_i_dot_vec(const double coef, const ElementField& vec, const CurrentFiniteElement& fePhi, const CurrentFiniteElement& fePsi, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 417 | |||
| 418 | //---------------------------------------------------------------------- | ||
| 419 | //! Operator for DG on faces $ \int_{\Sigma} eps_psi_j_dot_vec_phi_i $, | ||
| 420 | //! where $\Sigma$ is a boundary element. | ||
| 421 | //--------------------------------------------------------------------- | ||
| 422 | void eps_psi_j_dot_vec_phi_i(const double coef, const ElementField& vec, const CurrentFiniteElement& fePhi, const CurrentFiniteElement& fePsi, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 423 | |||
| 424 | //---------------------------------------------------------------------- | ||
| 425 | //! Operator for DG on faces $ \int_{\Sigma} grad_psi_j_dot_vec_phi_i $, | ||
| 426 | //! where $\Sigma$ is a boundary element. | ||
| 427 | //--------------------------------------------------------------------- | ||
| 428 | void grad_psi_j_dot_vec_phi_i(const double coef, const ElementField& vec, const CurrentFiniteElement& fePhi, const CurrentFiniteElement& fePsi, const std::size_t iblock, const std::size_t jblock, const std::size_t num); | ||
| 429 | |||
| 430 | //---------------------------------------------------------------------- | ||
| 431 | //! Operator for DG on faces $ \int_{\Sigma} f_dot_n_psi_j_phi_i $, | ||
| 432 | //! where $\Sigma$ is a boundary element. | ||
| 433 | // elt_f is the function f | ||
| 434 | // elt_n is the "normal". This is a constant elementField with as many components as elt_f | ||
| 435 | //--------------------------------------------------------------------- | ||
| 436 | void f_dot_n_psi_j_phi_i(const double coef, const CurrentFiniteElement& fePhi, const CurrentFiniteElement& fePsi, const ElementField& elt_f, const ElementField& elt_n, const std::size_t iblock, const std::size_t jblock, const std::size_t num, felInt stabFlag); | ||
| 437 | |||
| 438 | void const_phi_i_dot_n(const double coeff,const CurvilinearFiniteElement& fePhi, /*const CurvilinearFiniteElement& fePres,*/const std::size_t iBlock, const std::size_t jBlock/*, std::size_t num*/); | ||
| 439 | void const_phi_j_dot_n(const double coeff,const CurvilinearFiniteElement& fePhi, /*const CurvilinearFiniteElement& fePres,*/ const std::size_t iBlock, const std::size_t jBlock /*, std::size_t num*/); | ||
| 440 | |||
| 441 | void const_phi_i(const double coeff, const CurvilinearFiniteElement& fePhi, const std::size_t iBlock, const std::size_t jBlock); | ||
| 442 | void const_phi_j(const double coeff, const CurvilinearFiniteElement& fePhi, const std::size_t iBlock, const std::size_t jBlock); | ||
| 443 | |||
| 444 | |||
| 445 | |||
| 446 | ///@} | ||
| 447 | ///@name Access | ||
| 448 | ///@{ | ||
| 449 | |||
| 450 | 99951493 | UBlasMatrix& mat() { | |
| 451 | 99951493 | return m_mat; | |
| 452 | } | ||
| 453 | 8 | std::size_t numBlockRow() const { | |
| 454 | 8 | return m_numRow.size(); | |
| 455 | } | ||
| 456 | 8 | std::size_t numBlockCol() const { | |
| 457 | 8 | return m_numCol.size(); | |
| 458 | } | ||
| 459 | |||
| 460 | //std::size_t numRow() const {return m_numRow;} | ||
| 461 | |||
| 462 | //std::size_t numCol() const {return m_numCol;} | ||
| 463 | |||
| 464 | ///@} | ||
| 465 | ///@name Inquiry | ||
| 466 | ///@{ | ||
| 467 | |||
| 468 | ///@} | ||
| 469 | ///@name Input and output | ||
| 470 | ///@{ | ||
| 471 | |||
| 472 | void print(std::size_t verbose,std::ostream& c=std::cout); | ||
| 473 | |||
| 474 | ///@} | ||
| 475 | ///@name Friends | ||
| 476 | ///@{ | ||
| 477 | |||
| 478 | ///@} | ||
| 479 | protected: | ||
| 480 | ///@name Protected static Member Variables | ||
| 481 | ///@{ | ||
| 482 | |||
| 483 | ///@} | ||
| 484 | ///@name Protected member Variables | ||
| 485 | ///@{ | ||
| 486 | |||
| 487 | ///@} | ||
| 488 | ///@name Protected Operators | ||
| 489 | ///@{ | ||
| 490 | |||
| 491 | ///@} | ||
| 492 | ///@name Protected Operations | ||
| 493 | ///@{ | ||
| 494 | |||
| 495 | ///@} | ||
| 496 | ///@name Protected Access | ||
| 497 | ///@{ | ||
| 498 | |||
| 499 | ///@} | ||
| 500 | ///@name Protected Inquiry | ||
| 501 | ///@{ | ||
| 502 | |||
| 503 | ///@} | ||
| 504 | ///@name Protected LifeCycle | ||
| 505 | ///@{ | ||
| 506 | |||
| 507 | ///@} | ||
| 508 | private: | ||
| 509 | ///@name Private static Member Variables | ||
| 510 | ///@{ | ||
| 511 | |||
| 512 | ///@} | ||
| 513 | ///@name Private member Variables | ||
| 514 | ///@{ | ||
| 515 | |||
| 516 | UBlasMatrix m_mat; //!< The element matrix | ||
| 517 | std::vector<std::size_t> m_numRow; //!< m_numRow[i]=number of rows in the i-th block row | ||
| 518 | std::vector<std::size_t> m_firstRow; //!< m_firstRow[i]=index of the first row of i-th block row | ||
| 519 | std::vector<std::size_t> m_numCol; //!< m_numCol[i]=number of columns in the i-th block col | ||
| 520 | std::vector<std::size_t> m_firstCol; //!< m_firstCol[i]=index of the first column of i-th block col | ||
| 521 | |||
| 522 | // TODO:: All temp classes should be replaced with on time execution (otherwise much memory consumed) | ||
| 523 | UBlasMatrix m_tmpMat; //!< temporary element matrix (same size as m_mat) | ||
| 524 | UBlasVector m_tmpVecCoor;//!< temporary vector (size = numCoor) | ||
| 525 | UBlasVector m_tmpVec2Coor;//!< temporary vector (size = numCoor) | ||
| 526 | UBlasVector m_tmpVecDof;//!< temporary vector (size = numDof) | ||
| 527 | UBlasVector m_tmpVec2Dof;//!< temporary vector (size = numDof) | ||
| 528 | UBlasVector m_tmpVecQuadPoint; //!< temporary vector (size = numQuadraturePoint) | ||
| 529 | |||
| 530 | ///@} | ||
| 531 | ///@name Private Operators | ||
| 532 | ///@{ | ||
| 533 | |||
| 534 | ///@} | ||
| 535 | ///@name Private Operations | ||
| 536 | ///@{ | ||
| 537 | |||
| 538 | ///@} | ||
| 539 | ///@name Private Access | ||
| 540 | ///@{ | ||
| 541 | |||
| 542 | ///@} | ||
| 543 | ///@name Private Inquiry | ||
| 544 | ///@{ | ||
| 545 | |||
| 546 | ///@} | ||
| 547 | ///@name Private LifeCycle | ||
| 548 | ///@{ | ||
| 549 | |||
| 550 | ///@} | ||
| 551 | }; | ||
| 552 | ///@} | ||
| 553 | ///@name Type Definitions | ||
| 554 | ///@{ | ||
| 555 | |||
| 556 | ///@} | ||
| 557 | } /* namespace felisce.*/ | ||
| 558 | |||
| 559 | #endif | ||
| 560 |