| Directory: | ./ |
|---|---|
| File: | FiniteElement/elementField.cpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 145 | 333 | 43.5% |
| Branches: | 68 | 297 | 22.9% |
| 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-F. Gerbeau | ||
| 13 | // | ||
| 14 | |||
| 15 | // System includes | ||
| 16 | |||
| 17 | // External includes | ||
| 18 | |||
| 19 | // Project includes | ||
| 20 | #include "FiniteElement/elementField.hpp" | ||
| 21 | #include "PETScInterface/petscVector.hpp" | ||
| 22 | #include "PETScInterface/AbstractPETScObjectInterface.hpp" | ||
| 23 | |||
| 24 | namespace felisce | ||
| 25 | { | ||
| 26 | 218704 | ElementField::ElementField(): | |
| 27 | 218704 | m_numComp(0), | |
| 28 |
1/2✓ Branch 2 taken 218704 times.
✗ Branch 3 not taken.
|
218704 | m_dim(0) |
| 29 | { | ||
| 30 | 218704 | m_findDofPositionInPetscOrdering = nullptr; | |
| 31 | 218704 | } | |
| 32 | |||
| 33 | /***********************************************************************************/ | ||
| 34 | /***********************************************************************************/ | ||
| 35 | |||
| 36 | 218704 | ElementField::~ElementField() | |
| 37 | { | ||
| 38 |
2/2✓ Branch 0 taken 217539 times.
✓ Branch 1 taken 1165 times.
|
218704 | if (m_findDofPositionInPetscOrdering) { |
| 39 |
1/2✓ Branch 0 taken 217539 times.
✗ Branch 1 not taken.
|
217539 | delete [] m_findDofPositionInPetscOrdering; |
| 40 | 217539 | m_findDofPositionInPetscOrdering = nullptr; | |
| 41 | } | ||
| 42 | 218704 | } | |
| 43 | |||
| 44 | /***********************************************************************************/ | ||
| 45 | /***********************************************************************************/ | ||
| 46 | |||
| 47 | 6622 | void ElementField::initialize(const ElementFieldType type, const int numComp) | |
| 48 | { | ||
| 49 | 6622 | m_type = type; | |
| 50 | 6622 | m_numComp = numComp; | |
| 51 | 6622 | m_dim = 1; | |
| 52 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6622 times.
|
6622 | FEL_ASSERT(m_type == CONSTANT_FIELD); |
| 53 | 6622 | val.clear(); | |
| 54 | 6622 | val.resize(m_numComp,1); | |
| 55 | 6622 | val.clear(); | |
| 56 | 6622 | } | |
| 57 | |||
| 58 | /***********************************************************************************/ | ||
| 59 | /***********************************************************************************/ | ||
| 60 | |||
| 61 | 45188 | void ElementField::initialize(const ElementFieldType type, const CurrentFiniteElement& fe,const int numComp) | |
| 62 | { | ||
| 63 | 45188 | m_type = type; | |
| 64 | 45188 | m_numComp = numComp; | |
| 65 |
3/4✓ Branch 0 taken 12 times.
✓ Branch 1 taken 30270 times.
✓ Branch 2 taken 14906 times.
✗ Branch 3 not taken.
|
45188 | switch (m_type) { |
| 66 | 12 | case CONSTANT_FIELD: | |
| 67 | 12 | m_dim = 1; | |
| 68 | 12 | break; | |
| 69 | 30270 | case DOF_FIELD: | |
| 70 | 30270 | m_dim = fe.numDof(); | |
| 71 | 30270 | break; | |
| 72 | 14906 | case QUAD_POINT_FIELD: | |
| 73 | 14906 | m_dim = fe.numQuadraturePoint(); | |
| 74 | 14906 | break; | |
| 75 | // Default case should appear with a warning at compile time instead of an error in runtime | ||
| 76 | // (that's truly the point of using enums as switch cases) | ||
| 77 | // default: | ||
| 78 | //FEL_ERROR("Unknown element field type"); | ||
| 79 | // break; | ||
| 80 | } | ||
| 81 | 45188 | val.clear(); | |
| 82 | 45188 | val.resize(m_numComp,m_dim); | |
| 83 | 45188 | val.clear(); | |
| 84 | |||
| 85 |
2/2✓ Branch 0 taken 44308 times.
✓ Branch 1 taken 880 times.
|
45188 | if (m_findDofPositionInPetscOrdering) { |
| 86 |
1/2✓ Branch 0 taken 44308 times.
✗ Branch 1 not taken.
|
44308 | delete [] m_findDofPositionInPetscOrdering; |
| 87 | 44308 | m_findDofPositionInPetscOrdering = nullptr; | |
| 88 | } | ||
| 89 |
1/2✓ Branch 1 taken 45188 times.
✗ Branch 2 not taken.
|
45188 | m_findDofPositionInPetscOrdering = new felInt[fe.numDof()*m_numComp]; |
| 90 | 45188 | } | |
| 91 | |||
| 92 | /***********************************************************************************/ | ||
| 93 | /***********************************************************************************/ | ||
| 94 | |||
| 95 | 1239913 | void ElementField::initialize(const ElementFieldType type,const CurvilinearFiniteElement& fe,const int numComp) | |
| 96 | { | ||
| 97 | 1239913 | m_type = type; | |
| 98 | 1239913 | m_numComp = numComp; | |
| 99 |
3/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1235779 times.
✓ Branch 2 taken 4132 times.
✗ Branch 3 not taken.
|
1239913 | switch (m_type) { |
| 100 | 2 | case CONSTANT_FIELD: | |
| 101 | 2 | m_dim = 1; | |
| 102 | 2 | break; | |
| 103 | 1235779 | case DOF_FIELD: | |
| 104 | 1235779 | m_dim = fe.numDof(); | |
| 105 | 1235779 | break; | |
| 106 | 4132 | case QUAD_POINT_FIELD: | |
| 107 | 4132 | m_dim = fe.numQuadraturePoint(); | |
| 108 | 4132 | break; | |
| 109 | // Default case should appear with a warning at compile time instead of an error in runtime | ||
| 110 | // (that's truly the point of using enums as switch cases) | ||
| 111 | // default: | ||
| 112 | //FEL_ERROR("Unknown element field type"); | ||
| 113 | |||
| 114 | } | ||
| 115 | 1239913 | val.clear(); | |
| 116 | 1239913 | val.resize(m_numComp,m_dim); | |
| 117 | 1239913 | val.clear(); | |
| 118 | |||
| 119 |
2/2✓ Branch 0 taken 1023258 times.
✓ Branch 1 taken 216655 times.
|
1239913 | if (m_findDofPositionInPetscOrdering) { |
| 120 |
1/2✓ Branch 0 taken 1023258 times.
✗ Branch 1 not taken.
|
1023258 | delete [] m_findDofPositionInPetscOrdering; |
| 121 | 1023258 | m_findDofPositionInPetscOrdering = nullptr; | |
| 122 | } | ||
| 123 |
1/2✓ Branch 1 taken 1239913 times.
✗ Branch 2 not taken.
|
1239913 | m_findDofPositionInPetscOrdering = new felInt[fe.numDof()*m_numComp]; |
| 124 | 1239913 | } | |
| 125 | |||
| 126 | /***********************************************************************************/ | ||
| 127 | /***********************************************************************************/ | ||
| 128 | |||
| 129 | ✗ | void ElementField::initialize(const ElementFieldType type, const CurrentFiniteElementWithBd& fewbd,const int numComp) | |
| 130 | { | ||
| 131 | ✗ | m_type = type; | |
| 132 | ✗ | m_numComp = numComp; | |
| 133 | ✗ | switch (m_type) { | |
| 134 | ✗ | case CONSTANT_FIELD: | |
| 135 | ✗ | m_dim = 1; | |
| 136 | ✗ | break; | |
| 137 | ✗ | case DOF_FIELD: | |
| 138 | ✗ | m_dim = fewbd.numDof(); | |
| 139 | ✗ | break; | |
| 140 | ✗ | case QUAD_POINT_FIELD: | |
| 141 | ✗ | m_dim = fewbd.numQuadraturePointInternAndBd(); | |
| 142 | ✗ | break; | |
| 143 | // Default case should appear with a warning at compile time instead of an error in runtime | ||
| 144 | // (that's truly the point of using enums as switch cases) | ||
| 145 | // default: | ||
| 146 | // FEL_ERROR("Unknown element field type"); | ||
| 147 | } | ||
| 148 | ✗ | val.clear(); | |
| 149 | ✗ | val.resize(m_numComp,m_dim); | |
| 150 | ✗ | val.clear(); | |
| 151 | |||
| 152 | ✗ | if (m_findDofPositionInPetscOrdering) { | |
| 153 | ✗ | delete [] m_findDofPositionInPetscOrdering; | |
| 154 | ✗ | m_findDofPositionInPetscOrdering = nullptr; | |
| 155 | } | ||
| 156 | ✗ | m_findDofPositionInPetscOrdering = new felInt[fewbd.numDof()*m_numComp]; | |
| 157 | } | ||
| 158 | |||
| 159 | /***********************************************************************************/ | ||
| 160 | /***********************************************************************************/ | ||
| 161 | |||
| 162 | 720 | void ElementField::print(int verbose,std::ostream& c) | |
| 163 | { | ||
| 164 |
1/2✓ Branch 0 taken 720 times.
✗ Branch 1 not taken.
|
720 | if(verbose) { |
| 165 |
1/4✓ Branch 0 taken 720 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
720 | switch (m_type) { |
| 166 | 720 | case CONSTANT_FIELD: | |
| 167 | 720 | c << "Constant ElemField: "; | |
| 168 | 720 | break; | |
| 169 | ✗ | case DOF_FIELD: | |
| 170 | ✗ | c << "DOF field: " ; | |
| 171 | ✗ | break; | |
| 172 | ✗ | case QUAD_POINT_FIELD: | |
| 173 | ✗ | c << "Quadrature points ElemField: " ; | |
| 174 | ✗ | break; | |
| 175 | |||
| 176 | // Default case should appear with a warning at compile time instead of an error in runtime | ||
| 177 | // (that's truly the point of using enums as switch cases) | ||
| 178 | // default: | ||
| 179 | // FEL_ERROR("Unknown element field type"); | ||
| 180 | } | ||
| 181 | |||
| 182 |
1/2✓ Branch 0 taken 720 times.
✗ Branch 1 not taken.
|
720 | if(m_numComp == 1) { |
| 183 |
2/2✓ Branch 3 taken 720 times.
✓ Branch 4 taken 720 times.
|
1440 | for(int i=0; i<m_dim; i++) c << val(0,i) << " "; |
| 184 | } else { | ||
| 185 | ✗ | for(int icomp=0; icomp<m_numComp; icomp++) { | |
| 186 | ✗ | c << " Component " << icomp << ": "; | |
| 187 | ✗ | for(int i=0; i<m_dim; i++) c << val(icomp,i) << " "; | |
| 188 | } | ||
| 189 | } | ||
| 190 | 720 | c << std::endl; | |
| 191 | } | ||
| 192 | 720 | } | |
| 193 | |||
| 194 | /***********************************************************************************/ | ||
| 195 | /***********************************************************************************/ | ||
| 196 | |||
| 197 | 36161 | void ElementField::setValue(const std::vector<double>& value) | |
| 198 | { | ||
| 199 | // Now should do the same but also with quad_point_field. | ||
| 200 |
1/4✓ Branch 0 taken 36161 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
36161 | switch (m_type) { |
| 201 | 36161 | case CONSTANT_FIELD: | |
| 202 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 36161 times.
|
36161 | FEL_ASSERT(m_numComp == static_cast<int>(value.size())); |
| 203 |
2/2✓ Branch 0 taken 52989 times.
✓ Branch 1 taken 36161 times.
|
89150 | for (int i=0; i<m_numComp; i++) |
| 204 | 52989 | val(i,0) = value[i]; | |
| 205 | 36161 | break; | |
| 206 | ✗ | case QUAD_POINT_FIELD: | |
| 207 | ✗ | FEL_ASSERT(m_numComp == 1); | |
| 208 | ✗ | for (std::size_t i=0; i<value.size(); i++) | |
| 209 | ✗ | val(0,i) = value[i]; | |
| 210 | ✗ | break; | |
| 211 | ✗ | case DOF_FIELD: | |
| 212 | ✗ | FEL_ERROR("ElementField::setValue from std::vector is not yet implemented for this type of element field"); | |
| 213 | break; | ||
| 214 | } | ||
| 215 | 36161 | } | |
| 216 | |||
| 217 | /***********************************************************************************/ | ||
| 218 | /***********************************************************************************/ | ||
| 219 | |||
| 220 | 242 | void ElementField::setValue(const double value, const felInt icomp) | |
| 221 | { | ||
| 222 |
2/2✓ Branch 0 taken 726 times.
✓ Branch 1 taken 242 times.
|
968 | for (int k(0);k<m_dim;++k) { |
| 223 | 726 | val(icomp,k) = value; | |
| 224 | } | ||
| 225 | 242 | } | |
| 226 | |||
| 227 | /***********************************************************************************/ | ||
| 228 | /***********************************************************************************/ | ||
| 229 | |||
| 230 | 50262700 | void ElementField::setValue(const PetscVector& v,const CurBaseFiniteElement& fe, const felInt iel, | |
| 231 | const int idVariable, const AO& ao, const Dof& dof) | ||
| 232 | { | ||
| 233 | 50262700 | const int size = fe.numDof()*m_numComp; | |
| 234 |
1/4✓ Branch 0 taken 50262700 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
50262700 | switch (m_type) { |
| 235 | 50262700 | case DOF_FIELD: | |
| 236 | 50262700 | AOInterface::RetrieveDofPositionInPetscOrdering(m_findDofPositionInPetscOrdering, fe, iel, idVariable, ao, dof, m_numComp); | |
| 237 | 50262700 | v.getValues(size,m_findDofPositionInPetscOrdering,&val.data()[0]); | |
| 238 | 50262700 | break; | |
| 239 | |||
| 240 | ✗ | case QUAD_POINT_FIELD: { | |
| 241 | ✗ | felReal tmp[size]; | |
| 242 | ✗ | AOInterface::RetrieveDofPositionInPetscOrdering(m_findDofPositionInPetscOrdering, fe, iel, idVariable, ao, dof, m_numComp); | |
| 243 | ✗ | v.getValues(size,m_findDofPositionInPetscOrdering,tmp); | |
| 244 | |||
| 245 | ✗ | for(int iComp = 0; iComp<m_numComp; iComp++) { | |
| 246 | ✗ | for(int ig=0; ig<fe.numQuadraturePoint(); ig++) { | |
| 247 | ✗ | val(iComp,ig) = 0.; | |
| 248 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 249 | ✗ | val(iComp,ig) += tmp[iComp*fe.numDof() + iDof]*fe.phi[ig](iDof); | |
| 250 | } | ||
| 251 | } | ||
| 252 | } | ||
| 253 | ✗ | break; | |
| 254 | } | ||
| 255 | |||
| 256 | ✗ | case CONSTANT_FIELD: | |
| 257 | ✗ | FEL_ERROR("ElementField::setValue from Petsc std::vector is not yet implemented for this type of element field"); | |
| 258 | break; | ||
| 259 | } | ||
| 260 | 50262700 | } | |
| 261 | |||
| 262 | /***********************************************************************************/ | ||
| 263 | /***********************************************************************************/ | ||
| 264 | |||
| 265 | ✗ | void ElementField::setValue(const PetscVector& v,const CurvilinearFiniteElement& fe, const felInt iel, | |
| 266 | const int idVariable,ISLocalToGlobalMapping matisMapDof, const Dof& dof) | ||
| 267 | { | ||
| 268 | ✗ | const int size = fe.numDof()*m_numComp; | |
| 269 | ✗ | int* localIndexDof = new felInt[fe.numDof()*m_numComp]; | |
| 270 | ✗ | PetscInt output[size]; | |
| 271 | ✗ | int cpt = 0; | |
| 272 | felInt idDof; | ||
| 273 | ✗ | switch (m_type) { | |
| 274 | ✗ | case DOF_FIELD: | |
| 275 | ✗ | for(int iComp = 0; iComp<m_numComp; iComp++) { | |
| 276 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 277 | ✗ | dof.loc2glob(iel,iDof,idVariable,iComp,idDof); | |
| 278 | ✗ | localIndexDof[cpt] = idDof; | |
| 279 | ✗ | cpt++; | |
| 280 | } | ||
| 281 | } | ||
| 282 | ✗ | ISLocalToGlobalMappingApply(matisMapDof,size, localIndexDof, output); //m_findDofPositionInPetscOrdering); | |
| 283 | ✗ | v.getValues(size,output,&val.data()[0]); | |
| 284 | |||
| 285 | ✗ | break; | |
| 286 | ✗ | case QUAD_POINT_FIELD: { | |
| 287 | ✗ | FEL_ERROR("ElementField::setValue from Petsc std::vector is not yet implemented for this type of element field"); | |
| 288 | } | ||
| 289 | break; | ||
| 290 | ✗ | case CONSTANT_FIELD: | |
| 291 | ✗ | FEL_ERROR("ElementField::setValue from Petsc std::vector is not yet implemented for this type of element field"); | |
| 292 | break; | ||
| 293 | } | ||
| 294 | } | ||
| 295 | |||
| 296 | /***********************************************************************************/ | ||
| 297 | /***********************************************************************************/ | ||
| 298 | |||
| 299 | ✗ | void ElementField::setValueFdotNormal(const PetscVector& v,const CurvilinearFiniteElement& fe, const felInt iel, | |
| 300 | const int idVariable, const AO& ao, const Dof& dof) | ||
| 301 | { | ||
| 302 | //TODO: Attention, we assume a P0 normal!! | ||
| 303 | ✗ | FEL_ASSERT(m_numComp == 1); | |
| 304 | |||
| 305 | ✗ | const int size = fe.numDof()*fe.numCoor(); | |
| 306 | ✗ | std::vector<double> tmp(size); | |
| 307 | |||
| 308 | ✗ | felInt findDofPositionInPetscOrdering[size]; | |
| 309 | |||
| 310 | ✗ | int cpt = 0; | |
| 311 | felInt idDof; | ||
| 312 | ✗ | switch (m_type) { | |
| 313 | ✗ | case DOF_FIELD: | |
| 314 | ✗ | for(int iComp = 0; iComp<fe.numCoor(); iComp++) { | |
| 315 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 316 | ✗ | dof.loc2glob(iel,iDof,idVariable,iComp,idDof); | |
| 317 | ✗ | findDofPositionInPetscOrdering[cpt] = idDof; | |
| 318 | ✗ | cpt++; | |
| 319 | ✗ | val(0,iDof)=0; //I am not sure it is needed. | |
| 320 | } | ||
| 321 | } | ||
| 322 | ✗ | AOApplicationToPetsc(ao,size,findDofPositionInPetscOrdering); | |
| 323 | ✗ | v.getValues(size,findDofPositionInPetscOrdering,tmp.data()); | |
| 324 | ✗ | cpt = 0; | |
| 325 | ✗ | for(int iComp = 0; iComp<fe.numCoor(); iComp++) { | |
| 326 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 327 | ✗ | val(0,iDof) += tmp[cpt]*fe.normal[0](iComp); | |
| 328 | ✗ | cpt++; | |
| 329 | } | ||
| 330 | } | ||
| 331 | ✗ | break; | |
| 332 | ✗ | case QUAD_POINT_FIELD: | |
| 333 | case CONSTANT_FIELD: | ||
| 334 | ✗ | FEL_ERROR("ElementField::setValueFdotNormal from Petsc std::vector is not yet implemented for this type of element field"); | |
| 335 | break; | ||
| 336 | } | ||
| 337 | } | ||
| 338 | |||
| 339 | /***********************************************************************************/ | ||
| 340 | /***********************************************************************************/ | ||
| 341 | |||
| 342 | ✗ | void ElementField::setValue(const PetscVector& v,const CurvilinearFiniteElement& fe, const felInt iel, | |
| 343 | const int idVariable, const AO& ao, const Dof& dof, const int iComp) | ||
| 344 | { | ||
| 345 | ✗ | const int size = fe.numDof(); | |
| 346 | ✗ | switch (m_type) { | |
| 347 | ✗ | case DOF_FIELD: | |
| 348 | ✗ | AOInterface::RetrieveDofPositionInPetscOrderingSingleComponent(m_findDofPositionInPetscOrdering, fe, iel, idVariable, ao, dof, iComp); | |
| 349 | ✗ | v.getValues(size,m_findDofPositionInPetscOrdering,&val.data()[0]); | |
| 350 | ✗ | break; | |
| 351 | ✗ | case CONSTANT_FIELD: | |
| 352 | case QUAD_POINT_FIELD: | ||
| 353 | ✗ | FEL_ERROR("ElementField::setValue from Petsc std::vector is not yet implemented for this type of element field"); | |
| 354 | break; | ||
| 355 | } | ||
| 356 | } | ||
| 357 | |||
| 358 | /***********************************************************************************/ | ||
| 359 | /***********************************************************************************/ | ||
| 360 | |||
| 361 | ✗ | void ElementField::setValue(const PetscVector& v,const CurvilinearFiniteElement& fe, const felInt iel, | |
| 362 | const int idVariable, const AO& ao, const Dof& dof, const int iComp, const int efComp) | ||
| 363 | { | ||
| 364 | ✗ | const int size = fe.numDof(); | |
| 365 | ✗ | std::vector<double> tmp(size); | |
| 366 | ✗ | switch (m_type) { | |
| 367 | ✗ | case DOF_FIELD: | |
| 368 | ✗ | AOInterface::RetrieveDofPositionInPetscOrderingSingleComponent(m_findDofPositionInPetscOrdering, fe, iel, idVariable, ao, dof, iComp); | |
| 369 | ✗ | v.getValues(size,m_findDofPositionInPetscOrdering,&tmp[0]); | |
| 370 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 371 | ✗ | val(efComp,iDof) = tmp[iDof]; | |
| 372 | } | ||
| 373 | ✗ | break; | |
| 374 | ✗ | case CONSTANT_FIELD: | |
| 375 | case QUAD_POINT_FIELD: | ||
| 376 | ✗ | FEL_ERROR("ElementField::setValue from Petsc std::vector is not yet implemented for this type of element field"); | |
| 377 | break; | ||
| 378 | } | ||
| 379 | } | ||
| 380 | |||
| 381 | /***********************************************************************************/ | ||
| 382 | /***********************************************************************************/ | ||
| 383 | |||
| 384 | 405472 | void ElementField::setValue(const PetscVector& v,const CurrentFiniteElementWithBd& fe, felInt ibd, const felInt iel, | |
| 385 | const int idVariable, const AO& ao, const Dof& dof) | ||
| 386 | { | ||
| 387 |
1/2✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
|
405472 | CurvilinearFiniteElement* feCurv = fe.ptrBdEle(ibd); |
| 388 | 405472 | const int size = feCurv->numDof()*m_numComp; | |
| 389 | 405472 | int cpt = 0; | |
| 390 | felInt idDof; | ||
| 391 |
1/2✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
|
405472 | const std::vector<int>& idPointBdEle = fe.pointBdEle(ibd); |
| 392 | |||
| 393 |
1/4✓ Branch 0 taken 405472 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
405472 | switch (m_type) { |
| 394 | 405472 | case DOF_FIELD: | |
| 395 |
2/2✓ Branch 0 taken 810944 times.
✓ Branch 1 taken 405472 times.
|
1216416 | for(int iComp = 0; iComp<m_numComp; iComp++) { |
| 396 |
2/2✓ Branch 1 taken 1621888 times.
✓ Branch 2 taken 810944 times.
|
2432832 | for(int iDof=0; iDof<feCurv->numDof(); iDof++) { |
| 397 |
1/2✓ Branch 2 taken 1621888 times.
✗ Branch 3 not taken.
|
1621888 | dof.loc2glob(iel,idPointBdEle[iDof],idVariable,iComp,idDof); |
| 398 | 1621888 | m_findDofPositionInPetscOrdering[cpt] = idDof; | |
| 399 | 1621888 | cpt++; | |
| 400 | } | ||
| 401 | } | ||
| 402 |
1/2✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
|
405472 | AOApplicationToPetsc(ao,size,m_findDofPositionInPetscOrdering); |
| 403 |
2/4✓ Branch 2 taken 405472 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
|
405472 | v.getValues(size,m_findDofPositionInPetscOrdering,&val.data()[0]); |
| 404 | 405472 | break; | |
| 405 | |||
| 406 | ✗ | case QUAD_POINT_FIELD: { | |
| 407 | // size is small so no new needed | ||
| 408 | ✗ | felReal tmp[size]; | |
| 409 | |||
| 410 | ✗ | for(int iComp = 0; iComp<m_numComp; iComp++) { | |
| 411 | ✗ | for(int iDof=0; iDof<feCurv->numDof(); iDof++) { | |
| 412 | ✗ | dof.loc2glob(iel,idPointBdEle[iDof],idVariable,iComp,idDof); | |
| 413 | ✗ | m_findDofPositionInPetscOrdering[cpt] = idDof; | |
| 414 | ✗ | cpt++; | |
| 415 | } | ||
| 416 | } | ||
| 417 | ✗ | AOApplicationToPetsc(ao,size,m_findDofPositionInPetscOrdering); | |
| 418 | ✗ | v.getValues(size,m_findDofPositionInPetscOrdering,tmp); | |
| 419 | |||
| 420 | ✗ | for(int iComp = 0; iComp<m_numComp; iComp++) { | |
| 421 | ✗ | for(int ig=0; ig<feCurv->numQuadraturePoint(); ig++) { | |
| 422 | ✗ | val(iComp,ig) = 0.; | |
| 423 | ✗ | for(int iDof=0; iDof<feCurv->numDof(); iDof++) { | |
| 424 | ✗ | val(iComp,ig) += tmp[iComp*feCurv->numDof() + iDof]*feCurv->phi[ig](iDof); | |
| 425 | } | ||
| 426 | } | ||
| 427 | } | ||
| 428 | ✗ | break; | |
| 429 | } | ||
| 430 | |||
| 431 | ✗ | case CONSTANT_FIELD: | |
| 432 | ✗ | FEL_ERROR("ElementField::setValue from Petsc std::vector is not yet implemented for this type of element field"); | |
| 433 | break; | ||
| 434 | } | ||
| 435 | 405472 | } | |
| 436 | |||
| 437 | /***********************************************************************************/ | ||
| 438 | /***********************************************************************************/ | ||
| 439 | |||
| 440 | ✗ | void ElementField::setValue(const PetscVector& v,const CurBaseFiniteElement& fe, const felInt iel, | |
| 441 | const int idVariable, const Dof& dof) | ||
| 442 | { | ||
| 443 | felInt idLine; | ||
| 444 | double value; | ||
| 445 | felInt idDof; | ||
| 446 | ✗ | switch (m_type) { | |
| 447 | ✗ | case DOF_FIELD: | |
| 448 | ✗ | for(int iComp = 0; iComp<m_numComp; iComp++) { | |
| 449 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 450 | ✗ | dof.loc2glob(iel,iDof,idVariable,iComp,idDof); | |
| 451 | ✗ | idLine = idDof; | |
| 452 | ✗ | v.getValues(1,&idLine,&value); | |
| 453 | ✗ | val(iComp,iDof) = value; | |
| 454 | } | ||
| 455 | } | ||
| 456 | ✗ | break; | |
| 457 | ✗ | case CONSTANT_FIELD: | |
| 458 | case QUAD_POINT_FIELD: | ||
| 459 | ✗ | FEL_ERROR("ElementField::setValue from Petsc vector is not yet implemented for this type of element field"); | |
| 460 | break; | ||
| 461 | } | ||
| 462 | } | ||
| 463 | |||
| 464 | /***********************************************************************************/ | ||
| 465 | /***********************************************************************************/ | ||
| 466 | |||
| 467 | 8192 | void ElementField::setShape(CurrentFiniteElement fe, const ElementFieldType type, const int numComp) | |
| 468 | { | ||
| 469 | 8192 | m_type = type; | |
| 470 | 8192 | m_numComp = numComp; | |
| 471 | |||
| 472 |
1/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 8192 times.
✗ Branch 3 not taken.
|
8192 | switch(m_type) { |
| 473 | ✗ | case CONSTANT_FIELD: | |
| 474 | ✗ | m_dim = 1; | |
| 475 | ✗ | break; | |
| 476 | |||
| 477 | ✗ | case DOF_FIELD: | |
| 478 | ✗ | m_dim = fe.numDof(); | |
| 479 | ✗ | break; | |
| 480 | |||
| 481 | 8192 | case QUAD_POINT_FIELD: | |
| 482 | 8192 | m_dim = fe.numQuadraturePoint(); | |
| 483 | 8192 | break; | |
| 484 | // Default case should appear with a warning at compile time instead of an error in runtime | ||
| 485 | // (that's truly the point of using enums as switch cases) | ||
| 486 | // default: | ||
| 487 | // FEL_ERROR("Unexpected ElementField type.") | ||
| 488 | |||
| 489 | } | ||
| 490 | |||
| 491 | 8192 | val.clear(); | |
| 492 | 8192 | val.resize(m_numComp,m_dim); | |
| 493 | 8192 | val.clear(); | |
| 494 | |||
| 495 |
2/2✓ Branch 0 taken 8188 times.
✓ Branch 1 taken 4 times.
|
8192 | if (m_findDofPositionInPetscOrdering) { |
| 496 |
1/2✓ Branch 0 taken 8188 times.
✗ Branch 1 not taken.
|
8188 | delete [] m_findDofPositionInPetscOrdering; |
| 497 | 8188 | m_findDofPositionInPetscOrdering = nullptr; | |
| 498 | } | ||
| 499 |
1/2✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
|
8192 | m_findDofPositionInPetscOrdering = new felInt[fe.numDof()*m_numComp]; |
| 500 | 8192 | } | |
| 501 | |||
| 502 | /***********************************************************************************/ | ||
| 503 | /***********************************************************************************/ | ||
| 504 | |||
| 505 | 8192 | void ElementField::setValue(CurrentFiniteElement fe, ElementFieldDynamicValue &efv, const double time) | |
| 506 | { | ||
| 507 | Component comp; | ||
| 508 | const Point *refPointList; | ||
| 509 | Point geoPoint; | ||
| 510 | 8192 | double x = 0.; | |
| 511 | 8192 | double y = 0.; | |
| 512 | 8192 | double z = 0.; | |
| 513 | |||
| 514 | // Prepare shape (realloc val) | ||
| 515 |
4/8✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8192 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8192 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8192 times.
✗ Branch 11 not taken.
|
8192 | setShape(fe, efv.elementFieldType(), efv.numComp()); |
| 516 | |||
| 517 | // CONSTANT_FIELD | ||
| 518 |
2/4✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 8192 times.
|
8192 | if (efv.elementFieldType() == CONSTANT_FIELD) { |
| 519 | |||
| 520 | ✗ | for (int icomp=0 ; icomp<efv.numComp() ; icomp++) { | |
| 521 | ✗ | comp = efv.allComp[icomp]; | |
| 522 | ✗ | FEL_ASSERT( efv.typeValueOfElementField(comp)==FROM_CONSTANT ); | |
| 523 | ✗ | val(icomp,0) = efv.constant(comp); | |
| 524 | } | ||
| 525 | ✗ | return; | |
| 526 | } | ||
| 527 | |||
| 528 | // DOF_FIELD and QUAD_POINT_FIELD | ||
| 529 |
2/2✓ Branch 0 taken 24576 times.
✓ Branch 1 taken 8192 times.
|
32768 | for (int idim=0 ; idim<m_dim ; ++idim) { |
| 530 | |||
| 531 | /* get point coordinates */ | ||
| 532 |
2/6✓ Branch 1 taken 24576 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 24576 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
24576 | switch (efv.elementFieldType()) { |
| 533 | ✗ | case DOF_FIELD: | |
| 534 | ✗ | refPointList = fe.refEle().node(); | |
| 535 | ✗ | fe.coorMap(geoPoint, refPointList[idim]); | |
| 536 | ✗ | x = geoPoint.x(); | |
| 537 | ✗ | y = geoPoint.y(); | |
| 538 | ✗ | z = geoPoint.z(); | |
| 539 | ✗ | break; | |
| 540 | |||
| 541 | 24576 | case QUAD_POINT_FIELD: | |
| 542 |
2/2✓ Branch 1 taken 8192 times.
✓ Branch 2 taken 16384 times.
|
24576 | if ( !fe.hasQuadPtCoor() ) |
| 543 |
1/2✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
|
8192 | fe.computeCurrentQuadraturePoint(); |
| 544 | 24576 | x = fe.currentQuadPoint[idim].x(); | |
| 545 | 24576 | y = fe.currentQuadPoint[idim].y(); | |
| 546 | 24576 | z = fe.currentQuadPoint[idim].z(); | |
| 547 | 24576 | break; | |
| 548 | |||
| 549 | ✗ | case CONSTANT_FIELD: | |
| 550 | ✗ | FEL_ERROR("Wrong ElementFieldType value."); | |
| 551 | } | ||
| 552 | |||
| 553 | /* Set all component for this point */ | ||
| 554 |
3/4✓ Branch 1 taken 49152 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 24576 times.
✓ Branch 4 taken 24576 times.
|
49152 | for (int icomp=0 ; icomp<efv.numComp() ; icomp++) { |
| 555 | 24576 | comp = efv.allComp[icomp]; | |
| 556 | |||
| 557 |
2/6✓ Branch 1 taken 24576 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 24576 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
24576 | switch (efv.typeValueOfElementField(comp)) { |
| 558 | 24576 | case FROM_FUNCTION: | |
| 559 |
3/6✓ Branch 1 taken 24576 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24576 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24576 times.
✗ Branch 8 not taken.
|
24576 | val(icomp,idim) = efv.callbackXYZT(comp)(x,y,z,time); |
| 560 | 24576 | break; | |
| 561 | |||
| 562 | ✗ | case FROM_FILE: | |
| 563 | ✗ | FEL_ERROR("Not implemented."); | |
| 564 | break; | ||
| 565 | |||
| 566 | ✗ | case FROM_CONSTANT: | |
| 567 | ✗ | val(icomp,idim) = efv.constant(comp); | |
| 568 | ✗ | break; | |
| 569 | // Default case should appear with a warning at compile time instead of an error in runtime | ||
| 570 | // (that's truly the point of using enums as switch cases) | ||
| 571 | // default: | ||
| 572 | // FEL_ERROR("Wrong typeValueOfElementField value."); | ||
| 573 | } | ||
| 574 | } | ||
| 575 | } | ||
| 576 | } | ||
| 577 | |||
| 578 | /***********************************************************************************/ | ||
| 579 | /***********************************************************************************/ | ||
| 580 | |||
| 581 | ✗ | void ElementField::setValueMatching(double* DofValues,felInt* DofIdentities,felInt DofDimension, const CurvilinearFiniteElement& fe, const felInt iel, | |
| 582 | const int idVariable, const Dof& dof) | ||
| 583 | { | ||
| 584 | //DofValues = values of the variable to match | ||
| 585 | //DofIdentities = identity of the dofs in which to std::set the DofValues | ||
| 586 | ✗ | double value = 0.0; | |
| 587 | bool dofFound; | ||
| 588 | felInt idDof; | ||
| 589 | ✗ | switch (m_type) { | |
| 590 | ✗ | case DOF_FIELD: | |
| 591 | ✗ | for(int iDof=0; iDof<fe.numDof(); iDof++) { | |
| 592 | ✗ | for(int iComp = 0; iComp<m_numComp; iComp++) { | |
| 593 | ✗ | dof.loc2glob(iel,iDof,idVariable,iComp,idDof); | |
| 594 | ✗ | dofFound=false; | |
| 595 | ✗ | for(felInt i=0; i<DofDimension && dofFound == false; i++) { | |
| 596 | ✗ | if(DofIdentities[i]==idDof) { | |
| 597 | ✗ | value = DofValues[i]; | |
| 598 | ✗ | dofFound = true; | |
| 599 | } | ||
| 600 | } | ||
| 601 | ✗ | val(iComp,iDof) = value; | |
| 602 | } | ||
| 603 | } | ||
| 604 | ✗ | break; | |
| 605 | ✗ | case CONSTANT_FIELD: | |
| 606 | case QUAD_POINT_FIELD: | ||
| 607 | ✗ | FEL_ERROR("ElementField::setValueMatching is not yet implemented for this type of element field"); | |
| 608 | break; | ||
| 609 | } | ||
| 610 | } | ||
| 611 | |||
| 612 | /***********************************************************************************/ | ||
| 613 | /***********************************************************************************/ | ||
| 614 | |||
| 615 | 2377775 | UBlasMatrix& ElementField::get_val() | |
| 616 | { | ||
| 617 | 2377775 | return val; | |
| 618 | } | ||
| 619 | |||
| 620 | /***********************************************************************************/ | ||
| 621 | /***********************************************************************************/ | ||
| 622 | |||
| 623 | ✗ | UBlasVector ElementField::valAsVec() const | |
| 624 | { | ||
| 625 | ✗ | UBlasVector tmp(m_dim * m_numComp); | |
| 626 | ✗ | for ( int i(0);i<m_dim;++i) | |
| 627 | ✗ | for (int ic(0);ic<m_numComp;++ic) { | |
| 628 | ✗ | tmp(i + ic*m_numComp)=val(ic,i); | |
| 629 | } | ||
| 630 | ✗ | return tmp; | |
| 631 | } | ||
| 632 | |||
| 633 | } | ||
| 634 |