| Directory: | ./ |
|---|---|
| File: | Solver/linearProblem_template.hpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 47 | 128 | 36.7% |
| Branches: | 36 | 160 | 22.5% |
| 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 | ||
| 13 | // | ||
| 14 | |||
| 15 | #ifndef LINEARPROBLEM_TPP | ||
| 16 | #define LINEARPROBLEM_TPP | ||
| 17 | |||
| 18 | #ifndef LINEARPROBLEM_HPP | ||
| 19 | #error Include linearProblem.hpp instead | ||
| 20 | #endif | ||
| 21 | |||
| 22 | namespace felisce | ||
| 23 | { | ||
| 24 | /*! | ||
| 25 | \brief This function is a GENERIC ASSEMBLY LOOP on the boundary elements on a certain list of labels | ||
| 26 | |||
| 27 | \param[in] functionOfTheLoop is the function that actually do the computation element-wise | ||
| 28 | \param[in] labels a std::vector containing a list of integer labels where we want to loop | ||
| 29 | \param[in] a function to initialize for each element-type (typically to initialize the FE) | ||
| 30 | \param[in] updatefunction it could be implemented into the functionOfTheLoop, but like this is more general when dealing with duplicated supportDof | ||
| 31 | |||
| 32 | |||
| 33 | For several examples on how to use this class check linearProblemNSSimplifiedFSI.hpp | ||
| 34 | |||
| 35 | |||
| 36 | The template parameter theClass is going to be the class where the different functions are implemented. | ||
| 37 | |||
| 38 | // Nerd alert! | ||
| 39 | | I think that it is possible to avoid using a template by using some new binding techniques of c++11 to wrap | ||
| 40 | | the function of the class into functions not belonging to any particular class. (i.e. by a creating a free-function from a class function) | ||
| 41 | | | ||
| 42 | | The problem is that I can not simply declare void (*functionOfTheLoop)( felInt ), because the class of a | ||
| 43 | | a member function is like the first argument of the function itself and it therefore alters the signature (I guess..) | ||
| 44 | | | ||
| 45 | | The main drawback of this implementation is that I am not sure you can call this function directly with | ||
| 46 | | methods of different class, maybe if they belong to the same hierarchy is going to work. | ||
| 47 | // | ||
| 48 | */ | ||
| 49 | template<class theClass> | ||
| 50 | 1056 | void LinearProblem::assemblyLoopBoundaryGeneral( void (theClass::*functionOfTheLoop)( felInt ), | |
| 51 | const std::vector<felInt> & labels, | ||
| 52 | void (theClass::*initPerElementType)(), | ||
| 53 | void (theClass::*updateFunction)(const std::vector<Point*>&,const std::vector<int>&)) { | ||
| 54 | |||
| 55 | // some useful vector | ||
| 56 | 1056 | std::vector<Point*> elemPoint; | |
| 57 | 1056 | std::vector<felInt> elemIdPoint; | |
| 58 | 1056 | std::vector<felInt> vectorIdSupport; | |
| 59 | |||
| 60 | //to improve readability | ||
| 61 | 1056 | const std::vector<ElementType>& bagElementTypeDomainBoundary = m_meshLocal[m_currentMesh]->bagElementTypeDomainBoundary(); | |
| 62 | |||
| 63 | //initialization of the counter numElement | ||
| 64 | felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ]; //std::vector of size 23 number of the different type of elements | ||
| 65 |
2/2✓ Branch 0 taken 13200 times.
✓ Branch 1 taken 528 times.
|
27456 | for (int ityp=0; ityp < GeometricMeshRegion::m_numTypesOfElement; ityp++ ) { |
| 66 | 26400 | numElement[ityp]=0; | |
| 67 | } | ||
| 68 | |||
| 69 | //========loop on the element type of boundary (e.g. triangles..not a real loop) | ||
| 70 |
2/2✓ Branch 1 taken 528 times.
✓ Branch 2 taken 528 times.
|
2112 | for (std::size_t i = 0; i < bagElementTypeDomainBoundary.size(); ++i) { |
| 71 | |||
| 72 | 1056 | ElementType eltType = bagElementTypeDomainBoundary[i]; | |
| 73 |
1/2✓ Branch 1 taken 528 times.
✗ Branch 2 not taken.
|
1056 | defineCurvilinearFiniteElement(eltType); |
| 74 |
1/2✓ Branch 1 taken 528 times.
✗ Branch 2 not taken.
|
1056 | initElementArrayBD(); |
| 75 | |||
| 76 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 528 times.
✓ Branch 3 taken 528 times.
✗ Branch 4 not taken.
|
1056 | (static_cast<theClass*>(this)->*initPerElementType)(); |
| 77 | |||
| 78 | 1056 | int numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType]; | |
| 79 |
1/2✓ Branch 1 taken 528 times.
✗ Branch 2 not taken.
|
1056 | elemPoint.resize(numPointPerElt, nullptr); |
| 80 |
1/2✓ Branch 1 taken 528 times.
✗ Branch 2 not taken.
|
1056 | elemIdPoint.resize(numPointPerElt, 0); |
| 81 | |||
| 82 | //=======loop on subregion of the surface, with different labels. | ||
| 83 | 1056 | for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); | |
| 84 |
2/2✓ Branch 5 taken 1716 times.
✓ Branch 6 taken 528 times.
|
4488 | itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) { |
| 85 | 3432 | int currentLabel = itRef->first; | |
| 86 | 3432 | int numEltPerLabel = itRef->second.second; | |
| 87 | |||
| 88 |
3/4✓ Branch 4 taken 1716 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1068 times.
✓ Branch 8 taken 648 times.
|
3432 | if ( std::find( labels.begin(), labels.end(), currentLabel) != labels.end() ) { |
| 89 |
2/2✓ Branch 0 taken 654528 times.
✓ Branch 1 taken 1068 times.
|
1311192 | for ( felInt iel = 0; iel < numEltPerLabel; iel++) { |
| 90 | // return each id of point of the element and coordinate in two arrays: | ||
| 91 | // elemPoint and elemIdPoint. | ||
| 92 |
1/2✓ Branch 1 taken 654528 times.
✗ Branch 2 not taken.
|
1309056 | setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport); |
| 93 | |||
| 94 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 654528 times.
✓ Branch 3 taken 654528 times.
✗ Branch 4 not taken.
|
1309056 | (static_cast<theClass*>(this)->*updateFunction)(elemPoint,elemIdPoint); |
| 95 | |||
| 96 | // loop over all the support elements | ||
| 97 |
2/2✓ Branch 1 taken 654528 times.
✓ Branch 2 taken 654528 times.
|
2618112 | for (std::size_t it = 0; it < vectorIdSupport.size(); it++) { |
| 98 | // get the id of the support | ||
| 99 | 1309056 | felInt ielSupportDof = vectorIdSupport[it]; | |
| 100 | |||
| 101 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 654528 times.
✓ Branch 3 taken 654528 times.
✗ Branch 4 not taken.
|
1309056 | (static_cast<theClass*>(this)->*functionOfTheLoop)(ielSupportDof); |
| 102 | } | ||
| 103 | 1309056 | numElement[eltType]++; | |
| 104 | } | ||
| 105 | } else { | ||
| 106 | 1296 | numElement[eltType] += numEltPerLabel; | |
| 107 | } | ||
| 108 | } | ||
| 109 | } | ||
| 110 | } | ||
| 111 | |||
| 112 | |||
| 113 | /*! \brief this function assemble the mass matrix at the boundary. | ||
| 114 | * And the KSP | ||
| 115 | */ | ||
| 116 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 117 | template<class theClass> | ||
| 118 | ✗ | void LinearProblem::assembleMassBoundaryAndInitKSP( void (theClass::*functionOfTheLoop)( felInt ), | |
| 119 | const std::vector<felInt> & labels, | ||
| 120 | void (theClass::*initPerElementType)(), | ||
| 121 | void (theClass::*updateFunction)(const std::vector<Point*>&,const std::vector<int>&), | ||
| 122 | std::size_t iConn ) | ||
| 123 | { | ||
| 124 | ✗ | this->m_dofBD[iConn].allocateMatrixOnBoundary( m_massBD[iConn] ); | |
| 125 | ✗ | m_auxiliaryMatrix=m_massBD[iConn]; | |
| 126 | ✗ | m_auxiliaryDofBD = &(this->dofBD(iConn)); | |
| 127 | ✗ | if ( m_auxiliaryDofBD->hasDofsOnBoundary() ) { | |
| 128 | ✗ | this->assemblyLoopBoundaryGeneral( functionOfTheLoop, | |
| 129 | labels, | ||
| 130 | initPerElementType, | ||
| 131 | updateFunction); | ||
| 132 | ✗ | m_globPosColumn.clear(); | |
| 133 | ✗ | m_globPosRow.clear(); | |
| 134 | ✗ | m_matrixValues.clear(); | |
| 135 | |||
| 136 | ✗ | m_massBD[iConn].assembly(MAT_FINAL_ASSEMBLY); | |
| 137 | ✗ | if ( FelisceParam::instance().exportMassMatrix ) { | |
| 138 | ✗ | m_massBD[iConn].saveInBinaryFormat( this->m_dofBD[iConn].comm(), "Mass", FelisceParam::instance().resultDir); | |
| 139 | } | ||
| 140 | } | ||
| 141 | // Better to put it since m_auxiliaryMatrix will stay alive as long as the linear problem, | ||
| 142 | // while we would like to deallocate this matrix when the mass matrix gets out of scope | ||
| 143 | ✗ | m_auxiliaryMatrix.destroy(); | |
| 144 | ✗ | if ( m_auxiliaryDofBD->hasDofsOnBoundary() ) { | |
| 145 | ✗ | m_kspMassBD[iConn].init(this->dofBD(iConn).comm()); | |
| 146 | |||
| 147 | // Setting the type, the GMRES restart, the PC type, tolerances,... | ||
| 148 | ✗ | const auto& r_instance = FelisceParam::instance(this->instanceIndex()); | |
| 149 | ✗ | const std::string solver = r_instance.solver[m_identifier_solver]; | |
| 150 | ✗ | const std::string preconditioner = r_instance.preconditioner[m_identifier_solver]; | |
| 151 | ✗ | const double relativeTolerance = r_instance.relativeTolerance[m_identifier_solver]; | |
| 152 | ✗ | const double absoluteTolerance = r_instance.absoluteTolerance[m_identifier_solver]; | |
| 153 | ✗ | const int maxIteration = r_instance.maxIteration[m_identifier_solver]; | |
| 154 | ✗ | const int gmresRestart = r_instance.gmresRestart[m_identifier_solver]; | |
| 155 | ✗ | const bool initPrevSolution = r_instance.initSolverWithPreviousSolution[m_identifier_solver]; | |
| 156 | ✗ | const std::string preconOptions = r_instance.setPreconditionerOption[m_identifier_solver]; | |
| 157 | ✗ | m_kspMassBD[iConn].setKSPandPCType(solver, preconditioner); | |
| 158 | ✗ | m_kspMassBD[iConn].setTolerances(relativeTolerance, absoluteTolerance, maxIteration, gmresRestart); | |
| 159 | ✗ | m_kspMassBD[iConn].setKSPOptions(solver, initPrevSolution); | |
| 160 | ✗ | m_kspMassBD[iConn].setKSPOperator(m_massBD[iConn], preconOptions); | |
| 161 | } | ||
| 162 | } | ||
| 163 | |||
| 164 | |||
| 165 | /*! \brief this function assemble the mass matrix at the boundary. | ||
| 166 | * | ||
| 167 | */ | ||
| 168 | template<class theClass> | ||
| 169 | ✗ | void LinearProblem::assembleMassBoundaryOnly( void (theClass::*functionOfTheLoop)( felInt ), | |
| 170 | const std::vector<felInt> & labels, | ||
| 171 | void (theClass::*initPerElementType)(), | ||
| 172 | void (theClass::*updateFunction)(const std::vector<Point*>&,const std::vector<int>&), | ||
| 173 | DofBoundary &dofBD, PetscMatrix& massMatrix ) { | ||
| 174 | ✗ | dofBD.allocateMatrixOnBoundary( massMatrix ); | |
| 175 | ✗ | m_auxiliaryMatrix = massMatrix; | |
| 176 | ✗ | m_auxiliaryDofBD = &dofBD; | |
| 177 | ✗ | if ( dofBD.hasDofsOnBoundary() ) { | |
| 178 | ✗ | this->assemblyLoopBoundaryGeneral( functionOfTheLoop, | |
| 179 | labels, | ||
| 180 | initPerElementType, | ||
| 181 | updateFunction); | ||
| 182 | ✗ | m_globPosColumn.clear(); | |
| 183 | ✗ | m_globPosRow.clear(); | |
| 184 | ✗ | m_matrixValues.clear(); | |
| 185 | |||
| 186 | ✗ | massMatrix.assembly(MAT_FINAL_ASSEMBLY); | |
| 187 | } | ||
| 188 | // Better to put it since m_auxiliaryMatrix will stay alive as long as the linear problem, | ||
| 189 | // while we would like to deallocate this matrix when the mass matrix gets out of scope | ||
| 190 | ✗ | m_auxiliaryMatrix.destroy(); | |
| 191 | } | ||
| 192 | #endif | ||
| 193 | |||
| 194 | |||
| 195 | template <class FunctorXYZ> | ||
| 196 | inline void LinearProblem::set_U_0(const FunctorXYZ& functorXYZ, int iUnknown, int iComponent) { | ||
| 197 | felInt dof = 0; | ||
| 198 | Point pt; | ||
| 199 | std::set<felInt> allDof; | ||
| 200 | std::vector<double> value; | ||
| 201 | felInt sizeDof = 0; | ||
| 202 | felInt sizeUpdate = 0; | ||
| 203 | felInt iPosGlob; | ||
| 204 | double val; | ||
| 205 | |||
| 206 | int idVar = m_listUnknown.idVariable(iUnknown); | ||
| 207 | |||
| 208 | for ( felInt iel = 0; iel < (felInt)supportDofUnknown(iUnknown).iEle().size() - 1; iel++) { | ||
| 209 | for ( int iSupport = 0; iSupport < supportDofUnknown(iUnknown).getNumSupportDof(iel); iSupport++) { | ||
| 210 | m_dof.loc2glob(iel, iSupport, idVar, iComponent, dof); | ||
| 211 | allDof.insert(dof); | ||
| 212 | sizeUpdate = allDof.size(); | ||
| 213 | if ( sizeUpdate > sizeDof ) { | ||
| 214 | pt = supportDofUnknown(iUnknown).listNode()[supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[iel] + iSupport ]]; | ||
| 215 | val = functorXYZ(pt.x(), pt.y(), pt.z()); | ||
| 216 | iPosGlob = dof; | ||
| 217 | AOApplicationToPetsc(m_ao, 1, &iPosGlob); | ||
| 218 | m_seqSol.setValues( 1, &iPosGlob, &val, INSERT_VALUES); | ||
| 219 | m_sol.setValues( 1, &iPosGlob, &val, INSERT_VALUES); | ||
| 220 | } | ||
| 221 | sizeDof = allDof.size(); | ||
| 222 | } | ||
| 223 | } | ||
| 224 | |||
| 225 | m_seqSol.assembly(); | ||
| 226 | m_sol.assembly(); | ||
| 227 | } | ||
| 228 | |||
| 229 | template <class FunctorXYZ> | ||
| 230 | 8 | inline void LinearProblem::set_U_0_parallel(const FunctorXYZ& functorXYZ, int iUnknown, int iComponent) | |
| 231 | { | ||
| 232 | 8 | const int sizeLocal = m_numDofLocalUnknown[iUnknown]; | |
| 233 | 8 | const int idVar = m_listUnknown.idVariable(iUnknown); | |
| 234 | 8 | const int numCompOfVariable = m_listVariable[idVar].numComponent(); | |
| 235 | felInt indexGlobal; | ||
| 236 | Point pt; | ||
| 237 | |||
| 238 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
8 | std::vector<double> valueForPETSC(sizeLocal); |
| 239 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
8 | std::vector<felInt> idGlobalValue(sizeLocal,0); |
| 240 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
8 | std::vector<felInt> idLocalValue(sizeLocal); |
| 241 | 8 | std::iota(idLocalValue.begin(), idLocalValue.end(), 0); | |
| 242 | |||
| 243 |
1/2✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
|
8 | ISLocalToGlobalMappingApply(m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[iUnknown], sizeLocal, idLocalValue.data(), idGlobalValue.data()); |
| 244 | |||
| 245 |
2/2✓ Branch 0 taken 8450 times.
✓ Branch 1 taken 8 times.
|
8458 | for (felInt i = iComponent; i < sizeLocal; i+=numCompOfVariable) { |
| 246 | 8450 | indexGlobal = idGlobalValue[i]; | |
| 247 |
1/2✓ Branch 1 taken 8450 times.
✗ Branch 2 not taken.
|
8450 | AOPetscToApplication(m_ao,1,&indexGlobal); |
| 248 |
1/2✓ Branch 1 taken 8450 times.
✗ Branch 2 not taken.
|
8450 | findCoordinateWithIdDof(indexGlobal,pt); |
| 249 | 8450 | valueForPETSC[i] = functorXYZ(pt.x(), pt.y(), pt.z()); | |
| 250 | } | ||
| 251 | |||
| 252 |
1/2✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
8 | m_seqSol.setValues(sizeLocal, idGlobalValue.data(), valueForPETSC.data(), ADD_VALUES); |
| 253 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | m_seqSol.assembly(); |
| 254 |
1/2✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
8 | m_sol.setValues(sizeLocal, idGlobalValue.data(), valueForPETSC.data(), ADD_VALUES); |
| 255 |
1/2✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
|
8 | m_sol.assembly(); |
| 256 | 8 | } | |
| 257 | |||
| 258 | ✗ | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array) { | |
| 259 | Point pt; | ||
| 260 | ✗ | array.resize(m_numDof); | |
| 261 | felInt position; | ||
| 262 | ✗ | for( felInt i = 0; i < m_numDof; i++) { | |
| 263 | ✗ | findCoordinateWithIdDof(i,pt); | |
| 264 | ✗ | position = i; | |
| 265 | ✗ | AOApplicationToPetsc(m_ao,1,&position); | |
| 266 | ✗ | array[position] = func(pt.x(),pt.y(),pt.z()); | |
| 267 | } | ||
| 268 | } | ||
| 269 | |||
| 270 | // Warning: extraData size has to be numDof (not number of mesh points). | ||
| 271 | ✗ | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array, std::vector<double>& extraData) { | |
| 272 | Point pt; | ||
| 273 | ✗ | array.resize(m_numDof); | |
| 274 | felInt position; | ||
| 275 | ✗ | for( felInt i = 0; i < m_numDof; i++) { | |
| 276 | ✗ | findCoordinateWithIdDof(i,pt); | |
| 277 | ✗ | position = i; | |
| 278 | ✗ | AOApplicationToPetsc(m_ao,1,&position); | |
| 279 | ✗ | array[position] = func(pt.x(),pt.y(),pt.z(),extraData[i]); | |
| 280 | } | ||
| 281 | } | ||
| 282 | |||
| 283 | // Warning: extraData size has to be numDof (not number of mesh points). | ||
| 284 | ✗ | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func, std::vector<int>& array, std::vector<double>& extraData) { | |
| 285 | Point pt; | ||
| 286 | ✗ | array.resize(m_numDof); | |
| 287 | felInt position; | ||
| 288 | ✗ | for( felInt i = 0; i < m_numDof; i++) { | |
| 289 | ✗ | findCoordinateWithIdDof(i,pt); | |
| 290 | ✗ | position = i; | |
| 291 | ✗ | AOApplicationToPetsc(m_ao,1,&position); | |
| 292 | ✗ | array[position] = func(pt.x(),pt.y(),pt.z(),extraData[i]); | |
| 293 | } | ||
| 294 | } | ||
| 295 | |||
| 296 | ✗ | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func,double time, std::vector<double>& array) { | |
| 297 | Point pt; | ||
| 298 | ✗ | array.resize(m_numDof); | |
| 299 | felInt position; | ||
| 300 | ✗ | for( felInt i = 0; i < m_numDof; i++) { | |
| 301 | ✗ | findCoordinateWithIdDof(i,pt); | |
| 302 | ✗ | position = i; | |
| 303 | ✗ | AOApplicationToPetsc(m_ao,1,&position); | |
| 304 | ✗ | array[position] = func(pt.x(),pt.y(),pt.z(),time); | |
| 305 | } | ||
| 306 | } | ||
| 307 | |||
| 308 | // Warning: extraData size has to be numDof (not number of mesh points). | ||
| 309 | ✗ | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array, std::vector<double>& extraData,double time) { | |
| 310 | Point pt; | ||
| 311 | ✗ | array.resize(m_numDof); | |
| 312 | felInt position; | ||
| 313 | ✗ | for( felInt i = 0; i < m_numDof; i++) { | |
| 314 | ✗ | findCoordinateWithIdDof(i,pt); | |
| 315 | ✗ | position = i; | |
| 316 | ✗ | AOApplicationToPetsc(m_ao,1,&position); | |
| 317 | ✗ | array[position] = func(pt.x(),pt.y(),pt.z(),extraData[i],time); | |
| 318 | } | ||
| 319 | } | ||
| 320 | |||
| 321 | // Warning: extraData size has to be numDof (not number of mesh points). | ||
| 322 | ✗ | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func,double time, std::vector<double>& array, std::vector<double>& extraData) { | |
| 323 | Point pt; | ||
| 324 | ✗ | array.resize(m_numDof); | |
| 325 | felInt position; | ||
| 326 | ✗ | for( felInt i = 0; i < m_numDof; i++) { | |
| 327 | ✗ | findCoordinateWithIdDof(i,pt); | |
| 328 | ✗ | position = i; | |
| 329 | ✗ | AOApplicationToPetsc(m_ao,1,&position); | |
| 330 | ✗ | array[position] = func(pt.x(),pt.y(),pt.z(),time,extraData[i]); | |
| 331 | } | ||
| 332 | } | ||
| 333 | |||
| 334 | template <class Templ_functor> void LinearProblem::evalFunctionOnDof(const Templ_functor& func,int iComp, std::vector<double>& array) { | ||
| 335 | Point pt; | ||
| 336 | array.resize(m_numDof); | ||
| 337 | felInt position; | ||
| 338 | for( felInt i = 0; i < m_numDof; i++) { | ||
| 339 | findCoordinateWithIdDof(i,pt); | ||
| 340 | position = i; | ||
| 341 | AOApplicationToPetsc(m_ao,1,&position); | ||
| 342 | array[position] = func(iComp,pt.x(),pt.y(),pt.z()); | ||
| 343 | } | ||
| 344 | } | ||
| 345 | |||
| 346 | template <class Templ_functor> void LinearProblem::evalFunctionWithCompOnDof(const Templ_functor& func,std::vector<double>& array) { | ||
| 347 | Point pt; | ||
| 348 | for( felInt i = 0; i < m_numDof; i++) { | ||
| 349 | findCoordinateWithIdDof(i,pt); | ||
| 350 | for( int iComp = 0; iComp < 3; iComp++) { // change '3' to a m_numComp | ||
| 351 | array.push_back(func(iComp,pt.x(),pt.y(),pt.z())); | ||
| 352 | } | ||
| 353 | } | ||
| 354 | } | ||
| 355 | } | ||
| 356 | |||
| 357 | #endif // LINEAR_PROBLEM_TPP | ||
| 358 |