is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| Cal_SmallFemTermOfFemEquation.cpp (getdp-3.4.0-source.tgz) | : | Cal_SmallFemTermOfFemEquation.cpp (getdp-3.5.0-source.tgz) |
|---|
|
| // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege | | // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege |
| // | | // |
| // See the LICENSE.txt file for license information. Please report all | | // See the LICENSE.txt file for license information. Please report all |
| // issues on https://gitlab.onelab.info/getdp/getdp/issues. | | // issues on https://gitlab.onelab.info/getdp/getdp/issues. |
| // | | // |
| // Contributor(s): | | // Contributor(s): |
| // Nicolas Marsic | | // Nicolas Marsic |
| // | | // |
| | |
| #include "DofData.h" | | #include "DofData.h" |
| #include "Message.h" | | #include "Message.h" |
| | |
| skipping to change at line 25 | | skipping to change at line 25 |
|---|
| #include "SmallFem.h" | | #include "SmallFem.h" |
| #include "Mesh.h" | | #include "Mesh.h" |
| #include "FunctionSpace0Form.h" | | #include "FunctionSpace0Form.h" |
| #include "System.h" | | #include "System.h" |
| #include "SystemHelper.h" | | #include "SystemHelper.h" |
| #include "FormulationPoisson.h" | | #include "FormulationPoisson.h" |
| | |
| extern struct CurrentData Current; | | extern struct CurrentData Current; |
| | |
| // SmallFEM: helping functions // | | // SmallFEM: helping functions // |
|
| double fDirichlet0(fullVector<double>& xyz){ | | double fDirichlet0(fullVector<double> &xyz) { return -1; } |
| return -1; | | |
| } | | |
| | |
|
| double fDirichlet1(fullVector<double>& xyz){ | | double fDirichlet1(fullVector<double> &xyz) { return +1; } |
| return +1; | | |
| } | | |
| | |
|
| double fSource(fullVector<double>& xyz){ | | double fSource(fullVector<double> &xyz) { return 0; } |
| return 0; | | |
| } | | |
| | |
|
| void fMaterial(fullVector<double>& xyz, fullMatrix<double>& tensor){ | | void fMaterial(fullVector<double> &xyz, fullMatrix<double> &tensor) |
| | { |
| tensor.scale(0); | | tensor.scale(0); |
| | |
|
| if(xyz(0) > 0){ | | if(xyz(0) > 0) { |
| tensor(0, 0) = 1; | | tensor(0, 0) = 1; |
| tensor(1, 1) = 1; | | tensor(1, 1) = 1; |
| tensor(2, 2) = 1; | | tensor(2, 2) = 1; |
| } | | } |
|
| if(xyz(0) <= 0){ | | if(xyz(0) <= 0) { |
| tensor(0, 0) = 1; | | tensor(0, 0) = 1; |
| tensor(1, 1) = 1; | | tensor(1, 1) = 1; |
| tensor(2, 2) = 1; | | tensor(2, 2) = 1; |
| } | | } |
| } | | } |
| | |
|
| int getGoeIdFromGetDPId(GroupOfElement& goe, int getdpId){ | | int getGoeIdFromGetDPId(GroupOfElement &goe, int getdpId) |
| std::vector<const MElement*> element = goe.getAll(); | | { |
| | std::vector<const MElement *> element = goe.getAll(); |
| | |
| for(size_t i = 0; i < element.size(); i++) | | for(size_t i = 0; i < element.size(); i++) |
|
| if(element[i]->getNum() == getdpId) | | if(element[i]->getNum() == getdpId) return i; |
| return i; | | |
| | |
| throw Exception("Element not found!"); | | throw Exception("Element not found!"); |
| } | | } |
| | |
|
| struct Dof getGetDPDofFromSmallFEM(const sf::Dof& dofSF, | | struct Dof getGetDPDofFromSmallFEM(const sf::Dof &dofSF, |
| const DofManager<double>& dofM){ | | const DofManager<double> &dofM) |
| size_t globalId = dofM.getGlobalId(dofSF); | | { |
| | size_t globalId = dofM.getGlobalId(dofSF); |
| struct Dof dofDP; | | struct Dof dofDP; |
|
| dofDP.NumType = 1; | | dofDP.NumType = 1; |
| dofDP.Entity = dofSF.getEntity() + 1; | | dofDP.Entity = dofSF.getEntity() + 1; |
| dofDP.Harmonic = 0; | | dofDP.Harmonic = 0; |
| | |
|
| if(globalId == DofManager<double>::isFixedId()){ | | if(globalId == DofManager<double>::isFixedId()) { |
| dofDP.Type = 2; | | dofDP.Type = 2; |
| dofDP.Case.FixedAssociate.TimeFunctionIndex = 0; | | dofDP.Case.FixedAssociate.TimeFunctionIndex = 0; |
|
| dofDP.Val.s = dofM.getValue(dofSF); | | dofDP.Val.s = dofM.getValue(dofSF); |
| } | | } |
|
| else{ | | else { |
| dofDP.Type = 1; | | dofDP.Type = 1; |
| dofDP.Case.Unknown.NumDof = globalId + 1; | | dofDP.Case.Unknown.NumDof = globalId + 1; |
| dofDP.Case.Unknown.NonLocal = 0; | | dofDP.Case.Unknown.NonLocal = 0; |
| } | | } |
| | |
| return dofDP; | | return dofDP; |
| } | | } |
| | |
|
| void printDof(struct Dof* dof){ | | void printDof(struct Dof *dof) |
| std::cout << dof->NumType << ", " | | { |
| << dof->Entity << ", " | | std::cout << dof->NumType << ", " << dof->Entity << ", " << dof->Harmonic |
| << dof->Harmonic << ", " | | << ", " << dof->Type << "; "; |
| << dof->Type << "; "; | | |
| | | switch(dof->Type) { |
| switch(dof->Type){ | | case 1: |
| case 1: std::cout << dof->Case.Unknown.NumDof << ", " | | std::cout << dof->Case.Unknown.NumDof << ", " << dof->Case.Unknown.NonLocal; |
| << dof->Case.Unknown.NonLocal; break; | | break; |
| | | |
| case 2: std::cout << dof->Case.FixedAssociate.NumDof << ", " | | case 2: |
| << dof->Case.FixedAssociate.TimeFunctionIndex << ", " | | std::cout << dof->Case.FixedAssociate.NumDof << ", " |
| << dof->Val.s; break; | | << dof->Case.FixedAssociate.TimeFunctionIndex << ", " |
| | << dof->Val.s; |
| | break; |
| | |
| default: throw(Exception("Unknown GetDP Dof Type: %d", dof->Type)); | | default: throw(Exception("Unknown GetDP Dof Type: %d", dof->Type)); |
| } | | } |
| } | | } |
| | |
|
| void Cal_SmallFemTermOfFemEquation(struct Element* Element, | | void Cal_SmallFemTermOfFemEquation(struct Element *Element, |
| struct EquationTerm* EquationTerm_P, | | struct EquationTerm *EquationTerm_P, |
| struct QuantityStorage* QuantityStorage_P0){ | | struct QuantityStorage *QuantityStorage_P0) |
| extern int Flag_RHS; | | { |
| struct FemLocalTermActive* FI = EquationTerm_P->Case.LocalTerm.Active; | | extern int Flag_RHS; |
| | struct FemLocalTermActive *FI = EquationTerm_P->Case.LocalTerm.Active; |
| Current.flagAssDiag = 0; /*+++prov*/ | | Current.flagAssDiag = 0; /*+++prov*/ |
| | |
| /* treatment of MHBilinear-term in separate routine */ | | /* treatment of MHBilinear-term in separate routine */ |
|
| if(FI->MHBilinear){ | | if(FI->MHBilinear) { |
| /* if only the RHS of the system is to be calculated | | /* if only the RHS of the system is to be calculated |
| (in case of adaptive relaxation of the Newton-Raphson scheme) | | (in case of adaptive relaxation of the Newton-Raphson scheme) |
| the (expensive and redundant) calculation of this term can be skipped */ | | the (expensive and redundant) calculation of this term can be skipped */ |
|
| if(!Flag_RHS) | | if(!Flag_RHS) Message::Error("MHBilinear not in SmallFEM"); |
| Message::Error("MHBilinear not in SmallFEM"); | | |
| return; | | return; |
| } | | } |
| | |
| // Init SmallFEM (only once) // | | // Init SmallFEM (only once) // |
| static int once = 0; | | static int once = 0; |
|
| static Mesh* msh; | | static Mesh *msh; |
| static GroupOfElement* volume; | | static GroupOfElement *volume; |
| static GroupOfElement* boundary0; | | static GroupOfElement *boundary0; |
| static GroupOfElement* boundary1; | | static GroupOfElement *boundary1; |
| | | |
| static sf::FunctionSpace0Form* fs; | | static sf::FunctionSpace0Form *fs; |
| static FormulationPoisson* poisson; | | static FormulationPoisson *poisson; |
| static System<double>* sysPoisson; | | static System<double> *sysPoisson; |
| static const DofManager<double>* dofM; | | static const DofManager<double> *dofM; |
| | |
|
| static const std::vector<std::vector<sf::Dof> >* allDofField; | | static const std::vector<std::vector<sf::Dof> > *allDofField; |
| static const std::vector<std::vector<sf::Dof> >* allDofTest; | | static const std::vector<std::vector<sf::Dof> > *allDofTest; |
| | |
|
| if(!once){ | | if(!once) { |
| // Say it loudly and proudly! // | | // Say it loudly and proudly! // |
| Message::Direct("U s i n g S m a l l F E M . . ."); | | Message::Direct("U s i n g S m a l l F E M . . ."); |
| | |
| // Create a SmallFEM Formulation for Poisson // | | // Create a SmallFEM Formulation for Poisson // |
| // Get Domains // | | // Get Domains // |
|
| msh = new Mesh("mesh.msh"); | | msh = new Mesh("mesh.msh"); |
| volume = new GroupOfElement(msh->getFromPhysical(7)); | | volume = new GroupOfElement(msh->getFromPhysical(7)); |
| boundary0 = new GroupOfElement(msh->getFromPhysical(6)); | | boundary0 = new GroupOfElement(msh->getFromPhysical(6)); |
| boundary1 = new GroupOfElement(msh->getFromPhysical(5)); | | boundary1 = new GroupOfElement(msh->getFromPhysical(5)); |
| | |
| // Full Domain // | | // Full Domain // |
|
| std::vector<const GroupOfElement*> domain(3); | | std::vector<const GroupOfElement *> domain(3); |
| domain[0] = volume; | | domain[0] = volume; |
| domain[1] = boundary0; | | domain[1] = boundary0; |
| domain[2] = boundary1; | | domain[2] = boundary1; |
| | |
| // Get Order // | | // Get Order // |
| int order = 1; | | int order = 1; |
| | |
| // Function Space // | | // Function Space // |
| fs = new sf::FunctionSpace0Form(domain, order); | | fs = new sf::FunctionSpace0Form(domain, order); |
| | |
| | |
| skipping to change at line 170 | | skipping to change at line 168 |
|---|
| sysPoisson = new System<double>; | | sysPoisson = new System<double>; |
| sysPoisson->addFormulation(*poisson); | | sysPoisson->addFormulation(*poisson); |
| | |
| SystemHelper<double>::dirichlet(*sysPoisson, *fs, *boundary0, fDirichlet0); | | SystemHelper<double>::dirichlet(*sysPoisson, *fs, *boundary0, fDirichlet0); |
| SystemHelper<double>::dirichlet(*sysPoisson, *fs, *boundary1, fDirichlet1); | | SystemHelper<double>::dirichlet(*sysPoisson, *fs, *boundary1, fDirichlet1); |
| | |
| sysPoisson->assemble(); | | sysPoisson->assemble(); |
| dofM = &sysPoisson->getDofManager(); | | dofM = &sysPoisson->getDofManager(); |
| | |
| allDofField = &(poisson->field().getKeys(poisson->domain())); | | allDofField = &(poisson->field().getKeys(poisson->domain())); |
|
| allDofTest = &( poisson->test().getKeys(poisson->domain())); | | allDofTest = &(poisson->test().getKeys(poisson->domain())); |
| | |
| once = 1; | | once = 1; |
| } | | } |
| | |
| // Get GetDP Number of Dofs | | // Get GetDP Number of Dofs |
|
| //struct QuantityStorage* QuantityStorageEqu_P = FI->QuantityStorageEqu_P; | | // struct QuantityStorage* QuantityStorageEqu_P = FI->QuantityStorageEqu_P; |
| //struct QuantityStorage* QuantityStorageDof_P = FI->QuantityStorageDof_P; | | // struct QuantityStorage* QuantityStorageDof_P = FI->QuantityStorageDof_P; |
| //int Nbr_Dof = | | // int Nbr_Dof = |
| // FI->SymmetricalMatrix ? QuantityStorageEqu_P->NbrElementaryBasisFunction : | | // FI->SymmetricalMatrix ? QuantityStorageEqu_P->NbrElementaryBasisFunction : |
| // QuantityStorageDof_P->NbrElementaryBasisFunction; | | // QuantityStorageDof_P->NbrElementaryBasisFunction; |
|
| //int Nbr_Equ = | | // int Nbr_Equ = |
| // QuantityStorageEqu_P->NbrElementaryBasisFunction; | | // QuantityStorageEqu_P->NbrElementaryBasisFunction; |
| | |
| // SF Dofs | | // SF Dofs |
| int elementIdSF = getGoeIdFromGetDPId(*volume, Element->Num); | | int elementIdSF = getGoeIdFromGetDPId(*volume, Element->Num); |
| std::vector<sf::Dof> dofField = (*allDofField)[elementIdSF]; | | std::vector<sf::Dof> dofField = (*allDofField)[elementIdSF]; |
|
| std::vector<sf::Dof> dofTest = (*allDofTest)[elementIdSF]; | | std::vector<sf::Dof> dofTest = (*allDofTest)[elementIdSF]; |
| | |
| // Assemble | | // Assemble |
| int nDofField = dofField.size(); | | int nDofField = dofField.size(); |
|
| int nDofTest = dofTest.size(); | | int nDofTest = dofTest.size(); |
| | |
|
| for (int i = 0; i < nDofField; i++){ | | for(int i = 0; i < nDofField; i++) { |
| struct Dof dofI = getGetDPDofFromSmallFEM(dofField[i], *dofM); | | struct Dof dofI = getGetDPDofFromSmallFEM(dofField[i], *dofM); |
| | |
|
| for (int j = 0; j < nDofTest; j++){ | | for(int j = 0; j < nDofTest; j++) { |
| struct Dof dofJ = getGetDPDofFromSmallFEM(dofTest[j], *dofM); | | struct Dof dofJ = getGetDPDofFromSmallFEM(dofTest[j], *dofM); |
|
| double termIJ = poisson->weak(i, j, elementIdSF); | | double termIJ = poisson->weak(i, j, elementIdSF); |
| | |
|
| ((void (*)(struct Dof*, struct Dof*, double*)) | | ((void (*)(struct Dof *, struct Dof *, |
| FI->Function_AssembleTerm)(&dofI, &dofJ, &termIJ); | | double *))FI->Function_AssembleTerm)(&dofI, &dofJ, &termIJ); |
| } | | } |
| } | | } |
| | |
| // Done | | // Done |
| return; | | return; |
| } | | } |
| | |
| #else | | #else |
|
| void Cal_SmallFemTermOfFemEquation(struct Element* Element, | | void Cal_SmallFemTermOfFemEquation(struct Element *Element, |
| struct EquationTerm* EquationTerm_P, | | struct EquationTerm *EquationTerm_P, |
| struct QuantityStorage* QuantityStorage_P0){ | | struct QuantityStorage *QuantityStorage_P0) |
| | { |
| Message::Error("SmallFEM not activated"); | | Message::Error("SmallFEM not activated"); |
| } | | } |
| #endif | | #endif |
| | | |
| End of changes. 33 change blocks. |
|---|
| 79 lines changed or deleted | | 78 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 