diff --git a/applications/FSI/include/FSIassembly.hpp b/applications/FSI/include/FSIassembly.hpp index 282cc6fde..ae65a1ea3 100644 --- a/applications/FSI/include/FSIassembly.hpp +++ b/applications/FSI/include/FSIassembly.hpp @@ -163,7 +163,7 @@ void AssembleMatrixResFSI(MultiLevelProblem &ml_prob, unsigned level, const unsi vector indexVAR(2*dim+1); //vector indCOORD(dim); vector indVAR(3*dim+1); - vector SolType(3*dim+1); + vector SolType(3*dim+1); /// @todo is this really 3*d + 1 or 2*d+1 would be enough? for(unsigned ivar=0; ivar0: energy flows outside (cooling) QfluxDOTn<0: energy flows inside (heating) -void Temperature::heatflux_txyz(const double /*t*/, const double* /*xyz*/, double* qflux) const { - -// std::cout << "Temperature: Heatflux, check which coordinates are passed in here" << std::endl; -// Box* box= static_cast(_qtymap._phys._mesh->GetDomain()); -// const double thetaz = box->_domain_rtmap.get("thetaz"); - - qflux[0]=-2.1*0./**cos(thetaz)*/; - qflux[1]=0./**sin(thetaz)*/; - - return; - } - - - void Temperature::bc_flag_txyz(const double t, const double* xp, std::vector & bc_flag) const { // T' and its adjoint must be Dirichlet homogeneous everywhere on the boundary, by definition. diff --git a/applications/OptimalControl/fe_test/TempQuantities.hpp b/applications/OptimalControl/fe_test/TempQuantities.hpp index 46a528efa..561b4dd0b 100644 --- a/applications/OptimalControl/fe_test/TempQuantities.hpp +++ b/applications/OptimalControl/fe_test/TempQuantities.hpp @@ -20,10 +20,6 @@ class Temperature : public Quantity { void bc_flag_txyz(const double t, const double* xp, std::vector & flag) const; void initialize_xyz(const double* xp, std::vector & value) const; -//specific function - //this is the function of the IMPOSED DERIVATIVE of TEMPERATURE, aka heat flux - void heatflux_txyz(const double t,const double* xyz, double* qflux) const; - }; diff --git a/applications/OptimalControl/mhdopt/EqnMHD.cpp b/applications/OptimalControl/mhdopt/EqnMHD.cpp index c6f85721a..8659bd13b 100644 --- a/applications/OptimalControl/mhdopt/EqnMHD.cpp +++ b/applications/OptimalControl/mhdopt/EqnMHD.cpp @@ -71,7 +71,7 @@ void GenMatRhsMHD(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gr for (uint iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -155,7 +155,7 @@ void GenMatRhsMHD(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gr currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -359,7 +359,7 @@ for (uint fe = 0; fe < QL; fe++) { for (uint iel=0;iel < (nel_e - nel_b) ; iel++) { - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -449,7 +449,7 @@ for (uint fe = 0; fe < QL; fe++) { currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(ml_prob.GetMeshTwo()._iproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/mhdopt/EqnMHDAD.cpp b/applications/OptimalControl/mhdopt/EqnMHDAD.cpp index 308837d42..16c00ffae 100644 --- a/applications/OptimalControl/mhdopt/EqnMHDAD.cpp +++ b/applications/OptimalControl/mhdopt/EqnMHDAD.cpp @@ -67,7 +67,7 @@ const int NonStatMHDAD = (int) ml_prob.GetInputParser().get("NonStatMHDAD"); for (uint iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -140,7 +140,7 @@ const int NonStatMHDAD = (int) ml_prob.GetInputParser().get("NonStatMHDAD"); currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(ml_prob.GetMeshTwo()._iproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -321,7 +321,7 @@ for (uint fe = 0; fe < QL; fe++) { for (uint iel=0;iel < (nel_e - nel_b) ; iel++) { - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -364,7 +364,7 @@ for (uint fe = 0; fe < QL; fe++) { currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(ml_prob.GetMeshTwo()._iproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/mhdopt/EqnMHDCONT.cpp b/applications/OptimalControl/mhdopt/EqnMHDCONT.cpp index 3493736eb..609b375ae 100644 --- a/applications/OptimalControl/mhdopt/EqnMHDCONT.cpp +++ b/applications/OptimalControl/mhdopt/EqnMHDCONT.cpp @@ -77,7 +77,7 @@ using namespace femus; for (uint iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -138,7 +138,7 @@ using namespace femus; currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(ml_prob.GetMeshTwo()._iproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -373,7 +373,7 @@ for (uint fe = 0; fe < QL; fe++) { for (uint iel=0;iel < (nel_e - nel_b) ; iel++) { - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -410,7 +410,7 @@ for (uint fe = 0; fe < QL; fe++) { currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(ml_prob.GetMeshTwo()._iproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/mhdopt/EqnNS.cpp b/applications/OptimalControl/mhdopt/EqnNS.cpp index 8cf5db62b..bdf609622 100644 --- a/applications/OptimalControl/mhdopt/EqnNS.cpp +++ b/applications/OptimalControl/mhdopt/EqnNS.cpp @@ -159,7 +159,7 @@ const int NonStatNS = (int) ml_prob.GetInputParser().get("NonStatNS"); for (uint iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -250,7 +250,7 @@ const int NonStatNS = (int) ml_prob.GetInputParser().get("NonStatNS"); currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -547,7 +547,7 @@ for (uint fe = 0; fe < QL; fe++) { for (uint iel=0; iel < (nel_e - nel_b) ; iel++) { - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -601,7 +601,7 @@ for (uint fe = 0; fe < QL; fe++) { currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/mhdopt/EqnNSAD.cpp b/applications/OptimalControl/mhdopt/EqnNSAD.cpp index 3dbc7dd2d..e30b5d992 100644 --- a/applications/OptimalControl/mhdopt/EqnNSAD.cpp +++ b/applications/OptimalControl/mhdopt/EqnNSAD.cpp @@ -66,7 +66,7 @@ const int NonStatNSAD = (int) ml_prob.GetInputParser().get("NonStatNSAD"); for (uint iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -139,7 +139,7 @@ const int NonStatNSAD = (int) ml_prob.GetInputParser().get("NonStatNSAD"); currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -327,7 +327,7 @@ for (uint fe = 0; fe < QL; fe++) { for (uint iel=0;iel < (nel_e - nel_b) ; iel++) { - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -367,7 +367,7 @@ for (uint fe = 0; fe < QL; fe++) { currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/mhdopt/OptLoop.cpp b/applications/OptimalControl/mhdopt/OptLoop.cpp index 0e49ee6ab..1bf442904 100644 --- a/applications/OptimalControl/mhdopt/OptLoop.cpp +++ b/applications/OptimalControl/mhdopt/OptLoop.cpp @@ -511,15 +511,24 @@ double ComputeIntegral (const uint Level, const MultiLevelMeshTwo* mesh, const S const uint mesh_vb = VV; Mesh *mymsh = eqn->GetMLProb()._ml_msh->GetLevel(Level); - CurrentElem currelem(Level,VV,eqn,*mesh,eqn->GetMLProb().GetElemType()); - currelem.SetMesh(mymsh); - CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,eqn->GetMLProb().GetQrule(currelem.GetDim())); - - // processor index - const uint myproc = mesh->_iproc; + const unsigned myproc = mymsh->processor_id(); // geometry ----- const uint space_dim = mesh->get_dim(); + + + double integral=0.; + + +//parallel sum + const uint nel_e = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level+1]; + const uint nel_b = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level]; + + for (uint iel=0; iel < (nel_e - nel_b); iel++) { + + CurrentElem currelem(iel,myproc,Level,VV,eqn,*mesh,eqn->GetMLProb().GetElemType()); + currelem.SetMesh(mymsh); + CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,eqn->GetMLProb().GetQrule(currelem.GetDim())); //========= DOMAIN MAPPING CurrentQuantity xyz(currgp); @@ -546,18 +555,11 @@ double ComputeIntegral (const uint Level, const MultiLevelMeshTwo* mesh, const S VelDes._qtyptr = eqn->GetMLProb().GetQtyMap().GetQuantity("Qty_DesVelocity"); VelDes.VectWithQtyFillBasic(); VelDes.Allocate(); - - double integral=0.; const uint el_ngauss = eqn->GetMLProb().GetQrule(currelem.GetDim()).GetGaussPointsNumber(); - -//parallel sum - const uint nel_e = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level+1]; - const uint nel_b = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level]; - - for (uint iel=0; iel < (nel_e - nel_b); iel++) { - - currelem.SetDofobjConnCoords(mesh->_iproc,iel); + + + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/mhdopt/OptQuantities.cpp b/applications/OptimalControl/mhdopt/OptQuantities.cpp index fc6d81aa9..b8adc619f 100644 --- a/applications/OptimalControl/mhdopt/OptQuantities.cpp +++ b/applications/OptimalControl/mhdopt/OptQuantities.cpp @@ -219,32 +219,6 @@ void Velocity::Function_txyz(const double t,const double* xp, double* func) cons } -//============================================================= -void Velocity::strain_txyz(const double t, const double* xyz,double strain[][DIMENSION]) const { - -//here, tau is a tensor, so tau dot n is a vector which in general has a NORMAL and a TANGENTIAL component - - const double Lref = _qtymap.GetInputParser()->get("Lref"); - double ILref = 1./Lref; - const double lye = _qtymap.GetMeshTwo()->GetDomain()->_domain_rtmap.get("lye"); - const double x=xyz[0]; - const double y=xyz[1]; -#if DIMENSION==3 - const double z=xyz[2]; -#endif - - - strain[0][0] = 0.; //ux,x - strain[0][1] = strain[1][0] = 0. ;//0.5*(uy,x+ux,y) - strain[1][1] = 0.*(-(lye*ILref-y)); //uy,y -#if (DIMENSION==3) - strain[0][2] = strain[2][0] = 0. ; //0.5*(uz,x+ux,z) - strain[1][2] = strain[2][1] = 0. ; //0.5*(uy,z+uz,y) - strain[2][2] = 0. ; //uz,z -#endif - -return; -} //============================================================= @@ -366,26 +340,7 @@ void Temperature::Function_txyz(const double t, const double* xp,double* temp) c } -// ================================================= -void Temperature::heatflux_txyz(const double t, const double* xyz, double* qflux) const { - //the coordinates (x,y,z,t) of the VOLUME domain are NON-DIMENSIONAL - //and the function value must be nondimensional as well - //-----Nonhomogeneous Neumann------- - // Qflux = - k grad(T) by definition -// QfluxDOTn>0: energy flows outside (cooling) QfluxDOTn<0: energy flows inside (heating) - -std::cout << "Temperature: Heatflux, check which coordinates are passed in here" << std::endl; - Box* box = static_cast(_qtymap.GetMeshTwo()->GetDomain()); - const double thetaz = box->_domain_rtmap.get("thetaz"); - - qflux[0]=+700.*cos(thetaz); - qflux[1]=+700.*sin(thetaz); - #if (DIMENSION==3) - qflux[2]=0.; - #endif - return; - } diff --git a/applications/OptimalControl/mhdopt/OptQuantities.hpp b/applications/OptimalControl/mhdopt/OptQuantities.hpp index 2d0a293ec..3cb29d85d 100644 --- a/applications/OptimalControl/mhdopt/OptQuantities.hpp +++ b/applications/OptimalControl/mhdopt/OptQuantities.hpp @@ -141,15 +141,6 @@ class Velocity : public Quantity { void bc_flag_txyz(const double t, const double* xp, std::vector & flag) const; void initialize_xyz(const double* xp, std::vector & value) const; - //specific function - //this is the STRAIN DERIVATIVE of VELOCITY, so it must stay here - //from the physical and also mathematical point of view - //the shape funcs of the same order as v will then be used and so on -//multi-dimensional arrays must have bounds for all dimensions except the first - void strain_txyz(const double t, const double* xp,double strain[][DIMENSION]) const; //TODO remove this DIMENSION and use double pointers or std vector - - - }; class VelocityAdj : public Quantity { @@ -192,11 +183,6 @@ class Temperature : public Quantity { void Function_txyz(const double t, const double* xp,double* temp) const; -//specific function - //this is the function of the IMPOSED DERIVATIVE of TEMPERATURE, aka heat flux - void heatflux_txyz(const double t,const double* xyz, double* qflux) const; - - }; diff --git a/applications/OptimalControl/tempopt/EqnNS.cpp b/applications/OptimalControl/tempopt/EqnNS.cpp index 807bbc29d..21759687e 100644 --- a/applications/OptimalControl/tempopt/EqnNS.cpp +++ b/applications/OptimalControl/tempopt/EqnNS.cpp @@ -83,7 +83,7 @@ for (int iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -143,7 +143,7 @@ currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -364,154 +364,7 @@ for (uint fe = 0; fe < QL; fe++) { }//END VOLUME - - // ***************************************************************** - // ***************************************************************** - // ***************************************************************** - - {//BEGIN BOUNDARY // ***************************************************************** - - const uint mesh_vb = BB; - - const uint nel_e = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level+1]; - const uint nel_b = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level]; - - for (uint iel=0; iel < (nel_e - nel_b) ; iel++) { - - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); - currelem.SetMesh(mymsh); - CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); - - -//=========INTERNAL QUANTITIES (unknowns of the equation) ================== - CurrentQuantity VelOld(currgp); - VelOld._qtyptr = my_system.GetUnknownQuantitiesVector()[QTYZERO]; //an alternative cannot exist, because it is an Unknown of This Equation - VelOld.VectWithQtyFillBasic(); - VelOld.Allocate(); - - const uint qtyzero_ord = VelOld._FEord; - const uint qtyzero_ndof = VelOld._ndof; -// Velocity* vel_castqtyptr = static_cast(VelOld._qtyptr); //casting for quantity-specific functions - -//========= - CurrentQuantity pressOld(currgp); - pressOld._qtyptr = my_system.GetUnknownQuantitiesVector()[QTYONE]; - pressOld.VectWithQtyFillBasic(); - pressOld.Allocate(); - - const uint qtyone_ord = pressOld._FEord; - const uint qtyone_ndof = pressOld._ndof; - - //order - const uint qtyZeroToOne_DofOffset = VelOld._ndof*VelOld._dim; - -//========= END INTERNAL QUANTITIES (unknowns of the equation) ================= - -//=========EXTERNAL QUANTITIES (couplings) ===== - //========= //DOMAIN MAPPING - CurrentQuantity xyz(currgp); //domain - xyz._dim = space_dim; - xyz._FEord = MESH_MAPPING_FE; - xyz._ndof = currelem.GetElemType(xyz._FEord)->GetNDofs(); - xyz.Allocate(); - - //==================Quadratic domain, auxiliary, must be QUADRATIC!!! ========== - CurrentQuantity xyz_refbox(currgp); - xyz_refbox._dim = space_dim; - xyz_refbox._FEord = MESH_ORDER; - xyz_refbox._ndof = myel->GetElementFaceDofNumber(ZERO_ELEM,ZERO_FACE,BIQUADR_FE); - xyz_refbox.Allocate(); - -//=== auxiliary Operators at the boundary -// double strainU_g['dimension']['dimension']; -// double strainUtrDn_g['dimension']; - - currelem.Mat().zero(); - currelem.Rhs().zero(); - - currelem.SetDofobjConnCoords(myproc,iel); - currelem.SetMidpoint(); - - currelem.ConvertElemCoordsToMappingOrd(xyz); - currelem.TransformElemNodesToRef(ml_prob._ml_msh->GetDomain(),&xyz_refbox._val_dofs[0]); - currelem.SetElDofsBc(); - - VelOld.GetElemDofs(); - pressOld.GetElemDofs(); - -//============ BC ======= - int press_fl = currelem.Bc_ComputeElementBoundaryFlagsFromNodalFlagsForPressure(VelOld,pressOld); -//========END BC============ - -//============================================================== -//================== GAUSS LOOP (qp loop) ====================== -//============================================================== - const uint el_ngauss = ml_prob.GetQrule(currelem.GetDim()).GetGaussPointsNumber(); - - for (uint qp=0; qp< el_ngauss; qp++) { - -//======= "COMMON SHAPE PART"============================ - for (uint fe = 0; fe < QL; fe++) { - currgp.SetPhiElDofsFEVB_g (fe,qp); - currgp.SetDPhiDxezetaElDofsFEVB_g (fe,qp); -} - - const double det = dt*currgp.JacVectBB_g(xyz); - const double dtxJxW_g = det*ml_prob.GetQrule(currelem.GetDim()).GetGaussWeight(qp); -//=======end "COMMON SHAPE PART"=================================== - -//-------- pressure============== - //"predof" VS "post gauss method" - //pay attention to this fact: even if we are in the equation to which pressure is associated - //(pressure is an unknown here!), we might prefer using the FUNCTION instead, - //in order to get some values that DO NOT CHANGE with the solution - //so i'm using the FUNCTION of an UNKNOWN, it is not an external quantity - //clearly, this is ALTERNATIVE to the command - - - xyz_refbox.val_g(); - pressOld._qtyptr->Function_txyz(time,&xyz_refbox._val_g[0],&pressOld._val_g[0]); - - VelOld.val_g(); - -//============================================================== -//========= FILLING ELEMENT MAT/RHS (i loop) ==================== -//============================================================== - for (uint i=0; i_KK->add_matrix(currelem.Mat(),currelem.GetDofIndices()); - my_system._LinSolver[Level]->_RESC->add_vector(currelem.Rhs(),currelem.GetDofIndices()); - - - } - // end of BDRYelement loop - - - } - // END BOUNDARY ****************************** - my_system._LinSolver[Level]->_KK->close(); my_system._LinSolver[Level]->_RESC->close(); diff --git a/applications/OptimalControl/tempopt/EqnT.cpp b/applications/OptimalControl/tempopt/EqnT.cpp index 643b0b4b0..578dde5af 100644 --- a/applications/OptimalControl/tempopt/EqnT.cpp +++ b/applications/OptimalControl/tempopt/EqnT.cpp @@ -29,7 +29,6 @@ // application #include "TempQuantities.hpp" -#include "../../../src/equations/CurrentElem.hpp" // The question is: WHERE is the ORDER of the VARIABLES established? @@ -62,9 +61,6 @@ // For the derivatives, I think the point is: you must pick the REAL dphidx in the SAME ORDER as you pick the CORRESPONDING DOFS. // Now, my point is: on a given row, are you sure that the code picks the correct dphidx? -//TODO what happens for STANDARD OUTPUTS? can we do in such a way that EVERYTHING is printed TO FILE? -// we should REDIRECT TO THE *SAME* FILE ALL THE std output and std errors of ALL THE LIBRARIES! - //NOW, PAY ATTENTION: The "iel" written as "iel=0; iel < (nel_e - nel_b);" is used for PICKING the CONNECTIVITY from the ELEMENT CONNECTIVITY MAP! // But, the iel as DofObject Index must be given in the correct form! // So, I will distinguish iel into iel_mesh and iel_DofObj: @@ -103,10 +99,6 @@ void GenMatRhsT(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gri //==== AUXILIARY ============== std::vector dphijdx_g(space_dim); std::vector dphiidx_g(space_dim); - //-----Nonhomogeneous Neumann------ - // Qflux = - k grad(T) by definition -// QfluxDOTn>0: energy flows outside (cooling) QfluxDOTn<0: energy flows inside (heating) - std::vector Qflux_g(space_dim); my_system._LinSolver[Level]->_KK->zero(); my_system._LinSolver[Level]->_RESC->zero(); @@ -122,7 +114,7 @@ void GenMatRhsT(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gri for (uint iel=0; iel < (nel_e - nel_b); iel++) { - CurrentElem currelem(Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); + CurrentElem currelem(iel,myproc,Level,VV,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); @@ -179,19 +171,13 @@ void GenMatRhsT(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gri currelem.Mat().zero(); currelem.Rhs().zero(); - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); currelem.TransformElemNodesToRef(ml_prob._ml_msh->GetDomain(),&xyz_refbox._val_dofs[0]); -//MY EQUATION -//the elements are, for every level: -// 1)DOF INDICES -// 2)BC FLAGS -// 3)BC VALUES -// 1) and 2) are taken in a single vector, 3) are considered separately - + currelem.SetElDofsBc(); Tempold.GetElemDofs(); @@ -199,25 +185,13 @@ void GenMatRhsT(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gri TAdj.GetElemDofs(); -// =============== -// Now the point is this: there are several functions of space -// which are expressed with respect to a reference frame -//ok, now that the dofs are filled for xyz_refbox, I can use the el_average -//Well, the alternative is to consider the elem in the refbox as - //either a Vect or a CurrentElem ! - //I could consider it as another element, but only with the geometrical part! - xyz_refbox.SetElemAverage(); - -int domain_flag = ElFlagControl(xyz_refbox._el_average,ml_prob._ml_msh); + int domain_flag = ElFlagControl(xyz_refbox._el_average,ml_prob._ml_msh); //==================== //===== FILL the DOFS of the EXTERNAL QUANTITIES: you must assure that for every Vect the quantity is set correctly // for every Vect it must be clear if it belongs to an equation or not, and which equation it belongs to; // this is usually made clear by the related QUANTITY. -// Now the main thing to check is the difference between Vect WITH QUANTITY and Vect WITHOUT QUANTITY. - // it is better to avoid using GetElDofs if the Vect is internal, only if external - //Do not use GetElDofs if you want to pick an intermediate dof... if ( vel._eqnptr != NULL ) vel.GetElemDofs(); else vel._qtyptr->FunctionDof(vel,time,&xyz_refbox._val_dofs[0]); @@ -251,7 +225,7 @@ for (uint fe = 0; fe < QL; fe++) { TAdj.val_g(); vel.val_g(); Tdes.val_g(); - + // always remember to get the dofs for the variables you use! // The point is that you fill the dofs with different functions... // you should need a flag to check if the dofs have been correctly filled @@ -402,133 +376,6 @@ for (uint fe = 0; fe < QL; fe++) { }//END VOLUME - { //BEGIN BOUNDARY // ***************************************************************** - - const uint mesh_vb = BB; - - const uint nel_e = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level+1]; - const uint nel_b = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level]; - - for (uint iel=0;iel < (nel_e - nel_b) ; iel++) { - - - CurrentElem currelem(Level,BB,&my_system,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); - currelem.SetMesh(mymsh); - CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); - - -//=========INTERNAL QUANTITIES (unknowns of the equation) ========= - CurrentQuantity Tempold(currgp); - Tempold._qtyptr = my_system.GetUnknownQuantitiesVector()[0]; - Tempold.VectWithQtyFillBasic(); - Tempold.Allocate(); - -//==================================== - CurrentQuantity Tlift(currgp); - Tlift._qtyptr = ml_prob.GetQtyMap().GetQuantity("Qty_TempLift");//_UnknownQuantitiesVector[1]; - Tlift.VectWithQtyFillBasic(); - Tlift.Allocate(); - -//===================================== - CurrentQuantity TAdj(currgp); - TAdj._qtyptr = ml_prob.GetQtyMap().GetQuantity("Qty_TempAdj");//_UnknownQuantitiesVector[2]; - TAdj.VectWithQtyFillBasic(); - TAdj.Allocate(); - -//=========EXTERNAL QUANTITIES (couplings) ===== - //========= //DOMAIN MAPPING - CurrentQuantity xyz(currgp); //no quantity - xyz._dim = space_dim; - xyz._FEord = MESH_MAPPING_FE; - xyz._ndof = currelem.GetElemType(xyz._FEord)->GetNDofs(); - xyz.Allocate(); - - //==================Quadratic domain, auxiliary, must be QUADRATIC!!! ========== - CurrentQuantity xyz_refbox(currgp); //no quantity - xyz_refbox._dim = space_dim; - xyz_refbox._FEord = MESH_ORDER; - xyz_refbox._ndof = myel->GetElementFaceDofNumber(ZERO_ELEM,ZERO_FACE,BIQUADR_FE); - xyz_refbox.Allocate(); - -//===============Tdes===================== - CurrentQuantity Tdes(currgp); - Tdes._qtyptr = ml_prob.GetQtyMap().GetQuantity("Qty_TempDes"); - Tdes.VectWithQtyFillBasic(); - Tdes.Allocate(); - -// ========================================== -// ========================================== - - currelem.Mat().zero(); - currelem.Rhs().zero(); - - currelem.SetDofobjConnCoords(myproc,iel); - currelem.SetMidpoint(); - - currelem.ConvertElemCoordsToMappingOrd(xyz); - currelem.TransformElemNodesToRef(ml_prob._ml_msh->GetDomain(),&xyz_refbox._val_dofs[0]); - - currelem.SetElDofsBc(); - - Tempold.GetElemDofs(); - Tlift.GetElemDofs(); - TAdj.GetElemDofs(); - - //============ FLAGS ================ -//in order to do the flag here, since it is a "TRUE" NATURAL BOUNDARY CONDITION, -//it only suffices that SOME OF THE NODES ARE with bc=1, AT LEAST ONE -int el_Neum_flag=0; - uint Neum_sum=0; - for (uint i=0; i < Tempold._ndof; i++) Neum_sum += currelem.GetBCDofFlag()[i]; - for (uint i=0; i < Tempold._ndof; i++) Neum_sum += currelem.GetBCDofFlag()[i + Tempold._ndof]; - if ( Neum_sum == 2*Tempold._ndof ) { el_Neum_flag=1; } - -//==================================== - - const uint el_ngauss = ml_prob.GetQrule(currelem.GetDim()).GetGaussPointsNumber(); - - for (uint qp=0; qp< el_ngauss; qp++) { - -//======= "COMMON SHAPE PART"============================ - for (uint fe = 0; fe < QL; fe++) { - currgp.SetPhiElDofsFEVB_g (fe,qp); - currgp.SetDPhiDxezetaElDofsFEVB_g (fe,qp); - } - const double det = dt*currgp.JacVectBB_g(xyz); - const double dtxJxW_g = det * ml_prob.GetQrule(currelem.GetDim()).GetGaussWeight(qp); -//=======end "COMMON SHAPE PART"=================================== - - xyz.val_g(); - - static_cast(ml_prob.GetQtyMap().GetQuantity("Qty_Temperature"))->heatflux_txyz(time,&xyz._val_g[0],&Qflux_g[0]); - - Tempold.val_g(); //For the penalty Dirichlet //i need this for interpolating the old function at the gauss point - - double QfluxDn_g=Math::dot( &Qflux_g[0],currgp.get_normal_ptr(),space_dim); - - for (uint i=0; i_KK->add_matrix(currelem.Mat(),currelem.GetDofIndices()); - my_system._LinSolver[Level]->_RESC->add_vector(currelem.Rhs(),currelem.GetDofIndices()); - - } - // end of BDRYelement loop - - - }//END BOUNDARY - my_system._LinSolver[Level]->_KK->close(); my_system._LinSolver[Level]->_RESC->close(); @@ -542,232 +389,6 @@ int el_Neum_flag=0; } -//====================== - //TODO here I have to put the right offset back - //This function must receive a BOUNDARY ELEMENT according to the common "FEMuS OFFSET NUMERATION LevxSubd", - //and it must yield the VOLUME ELEMENT NUMBER again according to the common "FEMuS OFFSET NUMERATION LevxSubd", - //so that it has embedded the information of "BELONGING TO SOME GIVEN SUBDOMAIN" - //We start from a Boundary element belonging to some subdomain, and we want to get the CORRESPONDING VOLUME ELEMENT NUMBER, - //which is the volume element number IN THE "FEMuS OFFSET NUMERATION LevxSubd" NUMERATION, - //so that once you have it you already know to which processor you belong! - // Now this vol_iel that we have is already ok, because it is in the ABSOLUTE FEMuS VOLUME ELEMENT ORDERING - // Actually, in the mesh file generated by GenCase we never wrote so far the list of numbers of the elements, - // we only needed to write the offsets, - //then we reconstruct the numbering later, - //then we reassociate that numbering to the connectivity. - //It is like the connectivities were "orphan" of their respective numbers, but actually the connectivities - //were printed exactly in the new Femus ordering, - //(as a matter of fact, the FEMuS mesh file forgets about the libmesh ordering but only considers the femus one), - //so actually their order INTRINSICALLY GIVES the Femus ordering. - - //Now there is another story: the _el_bdry_to_vol returns an element number in the - // Femus ABSOLUTE Element numbering. - // But, the ABSOLUTE element numbering is used to pick the CONNECTIVITIES. - //To go from the ABSOLUTE to the RELATIVE Femus element numbering, - //you have to do like this: - // first go from the RELATIVE boundary to the ABSOLUTE boundary - //then from the absolute boundary you get the ABSOLUTE VOLUME element number - //finally, from the ABSOLUTE VOLUME YOU WANT TO GET THE DofObject - //TODO the thing that does not seem very nice to me is that - //there is an _el_bdry_to_vol map for every LEVEL, - //but the elements are NOT numbered by level, but according to - //the ABSOLUTE element numbering - - - //this iel_femus is ABSOLUTE, but PER LEVEL - // the following is ABSOLUTE, not even per level, COMPLETELY ABSOLUTE - //maybe it would be better to have it absolute but per level -//====================== - - -//====================== -// Now the dof indices on the boundary must match with the -//dofs on the volume -//If with constant elements you do not have any dof on the boundary, -//then you have to think like this: -//What is the meaning of ENFORCING a BOUNDARY CONDITION VALUE? -//The meaning is that you enforce the TRACE of the VOLUME SHAPE FUNCTIONS. -//Then, of course, enforcing the TRACE means enforcing the VOLUME DOF, for constant FE. -//With nodes on the boundary, you are enforcing the DOF ITSELF, because -//among all the nodes there are some which are ONLY BOUNDARY NODES. -//Again, enforcing those boundary nodes can be seen as a form of ENFORCING the TRACE of the -//VOLUME SHAPE FUNCTIONS associated to those BOUNDARY NODES. -//In fact, you enforce the DOF values so that you "stretch" your TRACE at the boundary. - -//Now, the point here is that we have a VOLUME DOF, iel, associated to the CONSTANT FE, -//which may not belong to the current row actually -//Wait, so, do we need to have the VOLUME iel in the list of BOUNDARY DOF INDICES? -//"Boundary dof indices" are the -// "dof_indices that are involved in considering the current BOUNDARY ELEMENT" -//DOF_INDICES means "ROWS and COLUMNS" of the matrix - -//I dont care if i am in a BOUNDARY INTEGRAL and the dofs I involve -//are not "boundary dofs" strictly speaking... -//I am only interested in INVOLVING ALL THE DOFS THAT ARE INVOLVED . -//Now, I need a map that for every boundary element -//gives me the corresponding VOLUME element in the mesh. -//For every volume element, i can give you the children, which can be one or two... -//for every boundary element, there is only one father, -//so let us do a simple map that goes from CHILDREN to FATHER. -//That is not gonna be in terms of DOF but in terms of MESH. -//So we go from BOUNDARY iel to VOLUME iel. -//We need to do this in the Gencase because it is there -//where we can use the Libmesh functions - - - - - - - - //on the other hand consider that the penalty term does not have to appear whenever we are inside the domain -//the problem is that being inside the domain doesnt mean being Dirichlet, but here bc=0 may mean either being inside the domain -//or being a Dirichlet node. We want to put the penalty term where bc=0,but not inside the domain!!! -//so it takes a flag for setting dirichlet -//a flag for setting that we are inside the domain...no,but we already know that we are on the boundary, because -//the problem is that some nodes are BOTH VOLUME AND BOUNDARY. Or better, the SOME VOLUME NODES are also BOUNDARY NODES... -//no no - - - //this should be removed I think: yes it should, - //I mean its useless. With the PENALTY APPROACH you dont have a flag 0-1, but 0-\infty! - //therefore, you only have to distinguish two types of elements: - // - the ELEMENTS in which penalty=\infty, where the Dirichlet integral dominates over everything - // - the ELEMENTS with penalty=0, for which we know that bc=1 so the bc_alldofs in front is not needed! - //in the first ones the Dirichlet integral, AND ONLY IT, - //in the second ones the Neumann integral, TOGETHER WITH THE REST OF THE EQUATION - //instead in the NODAL DIRICHLET BCs, Neumann is still ELEMENTAL, so we must implement it - //with Neum_flag, NOT WITH bc_alldofs again! - //in conclusion, bc_alldofs must be eliminated FOREVER, because it's NODAL, and no nodal thing must be used for the boundary INTEGRALS! - //what remains is Neum_flag when you dont use penalty - //Well,actually we HAVE to use bc_alldofs[i] in case of NODAL DIRICHLET. In fact, in that case - //ALL THE ROWS OF THE LINEAR SYSTEM must be multiplied by bc_alldofs, because - // the bc_qldof flag means "in this row we are putting or not a Dirichlet nodal value" - // instead, el_Neum_flag means "For the ELEMENT we consider we must compute the integral" -// (and this computation will involve more than one row, so we cannot isolate the rows) -//clearly, this situation creates a little mismatch, because it may happen -//that the nodal bc-s are 1-1-0, -//but the integral should be computed over ALL the element -//it seems like you have to interpolate with three functions but you stop the sum to 2 instead. -//The missing term should appear in equation "i" in some column "j", but equation "i" in the nodal case -//has all zero except 1 on the diagonal -//Removing that node corresponds to TRUNCATING the INTERPOLATION of the SHAPE FUNCTIONS (i index) at that gauss point. -//That doesnt give any problem because the "coefficient that doesn't multiply anything" would have multiplied -// a zero, so no fear! - - - - - - - //PAY ATTENTION here because there may be a slight problem if some errors in the update of Tempold._val_g are there, -//like it happened with p_old for pressure in NS, which suggested us to use a FUNCTION instead of the previous value - -//el_Neum_flag and el_penalty are clearly ALTERNATIVE: -// el_Neum_flag == 1 el_penalty == 0 -// el_Neum_flag == 0 el_penalty == 10e20 -//This means that you shouldnt need multiplying by el_Neum_flag in the Neumann integral. -//Well actually I would like to keep some symmetry... -//In practice, the basis is writing the general equation... Then the issue is enforcing Dirichlet. -//In the NODAL case it is bc_alldofs that rules. -//In the PENALTY case it is el_penalty that rules - //So, in both cases computing el_Neum_flag is useless - //in the nodal case, el_Neum_flag is 1 when bc is 1 -//You may want to use el_Neum_flag to avoid adding the small Neum integral to the 10e20 row, -//if you're very fond of symmetry... -//yes, it is a minor thing but it is harmless -//well, actually it might be also good when the numbers of the heat flux are very high... -//So we'll leave it, it looks as stabilyzing... - -// T' should be equal to -T_0 -//The point is that T' is of course homogeneous at the boundary, by definition. -//The part that takes into account the boundary conditions is only T_0: -//now try to set the bc free for T': in this case you get exactly T' = -T_0; -// but this is not what you should get. -//Of course, if you have different boundary conditions for T' and T_0, -//they can never be equal... -//With the magnetic field we didnt have b' = - B_e -//now try to solve for different T_0 (change the equation, put the laplacian, get sthg else...) -//is it true that different T_0 of the SAME boundary conditions lead to the same solution T' + T_0 ? -//Given different lifting functions of the same boudnarty conditions, -//what is the condition sucj that we do not depend on the particular lifting function? -// For the magnetic field we had div B = 0, and we claimed that in that case -//the final sum was the same -//if we had div B != 0, the final result in terms of u and p was different -//maybe that case was particular only for Hartmann, which is a linear, simplified, MHD problem? - -//our previous claim was: the NS + MHD operators where in such a form that their action -// on various B_e was leading to results that were only dependent on the boundary conditions. -// We should check this thing more THEORETICALLY: considering how the operators act on the -//decomposition b+ B_e , where B_e is or not divergence-free. - -//in order to find the equivalence, i think that we have to check the OPERATOR we have. -//if the operator is not sensible to some class of functions, or better if the action of that -//operator does not "filter out" the particular choice of the lifting function -// and only retain the boundary condition values, -//then that is not a good operator for us. -// we need a sort of "independent operator", or "operator of boundary conditions" -//LAPLACIAN -//LAPLACIAN + TRNSPORT - -// The point is: if the solution of the operator for T is unique, -// then whatever T_0 you set you get another T', but the sum (T' + T_0) is unique. -// And, of course, if the boundary conditions are the same, the sum (T' + T_0) is the same. - -//Now with the lifting of B you ALSO needed that the lifting function was DIVERGENCE-FREE: -// otherwise, you were getting different solutions for (u,p). -//But, in this case you dont have any particular constraint on T_0 - -//For the temperature, you may say that if the solution of the state equations is unique -//for a given boundary condition, then whatever T_0 you get, the sum T' + T_0 is unique. - -// Numerically speaking, the only point would be to be sure that the matrix is diagonal dominant -// (coming from coercivity) for every lifting function: clearly we know that coercivity -// is fulfilled only for certain lifting functions, for a certain epsilon value wrt the viscosity. - -// PROBLEMA: Boundary conditions for the STATE, ADJOINT AND CONTROL variables -// Secondo la mia teoria, sia per T' sia per la sua aggiunta le condizioni al boundary devono essere FISSE. - -// Come e' giusto per un controllo d boundary, l'effetto prevalente deve avvenire vicino al boundary. -// Quindi il controllo di boundary agisce sul boundary del dominio, pertanto non puo' spingere molto all'interno. - -//$$$$$$$$$$$$$$$$ Come posso aumentare la spinta all'interno? - -//$$$$$$$$$$$$$ Il ruolo delle boundary conditions: slego T_0 al boundary, pero' le altre non sono slegate -// al boundary. Non slegherei T', magari proverei a slegare l'aggiunta, -// che da la spinta al controllo, che da la spinta allo stato - -//$$$$$$$$$$$$$$$ Provare a cambiare con il target, o con la condizione iniziale di Stato/Controllo - -//$$$$$$$$$$ E' chiaro che il controllo di boundary non e' cosi' efficace come il controllo distribuito - -//CONDIZIONI di NEUMANN con il metodo del lifting: come si fanno a mettere? -//Ad esempio nella parte non controllata vorrei poter mettere dei pezzi adiabatici. -//Allora qual e' il punto: -// l'imposizione del controllo e l'imposizione della condizione di neumann -// avverrebbero nello stesso modo, -// cioe' non fissando il valore in quel nodo. -// Allora come fa il sistema a capire se un pezzo di boundary della lifting function -// e' un pezzo di CONTROLLO o un pezzo di NEUMANN OMOGENEO ? - -// Aspetta: se tu nell'equazione della temperatura stai lasciando tutto libero e -// non stai fissando l'integrale di boundary, allora stai fissando tutto uguale a zero... -//Poi pero' il pezzo dell'integrale di boundary per T_0 va a finire nell'equazione del controllo -// quando fai la variazione delta T_0, quindi nell'equazione del controllo avresti -// l'integrale di boundary dell'aggiunta... -// e allora la domanda e': Qual e' la condizione di Neumann dell'aggiunta? - - - -//Poi, come si fa a tenere la temperatura positiva? Un algoritmo serio dovrebbe evitare temperature negative... -// Se pero' l'equazione e' adimensionale, i valori adimensionali possono venire anche negativi, -// dipende da come hai adimensionalizzato: l'importante e' che poi i valori finali siano positivi -// SE IO AGGIUNGO AI VINCOLI ANCHE UN VINCOLO DI DISUGUAGLIANZA DICENDO CHE LA TEMPERATURA SIA POSITIVA? - - - -// Succede questo (quasi paradossale): se la regione di controllo e' vicino al boundary control, -// converge fino ad alpha massimo molto piccolo. -// Se la regione di controllo e' lontana dal boundary control, allora converge fino ad un alpha max molto piu' grande. + #endif \ No newline at end of file diff --git a/applications/OptimalControl/tempopt/OptLoop.cpp b/applications/OptimalControl/tempopt/OptLoop.cpp index 76e8d41a1..71bfac0a9 100644 --- a/applications/OptimalControl/tempopt/OptLoop.cpp +++ b/applications/OptimalControl/tempopt/OptLoop.cpp @@ -89,7 +89,19 @@ double ComputeIntegral (const uint Level, const MultiLevelMeshTwo* mesh, const S const uint mesh_vb = VV; - CurrentElem currelem(Level,VV,eqn,*mesh,eqn->GetMLProb().GetElemType()); + + + + double integral = 0.; + + +//parallel sum + const uint nel_e = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level+1]; + const uint nel_b = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level]; + + for (uint iel=0; iel < (nel_e - nel_b); iel++) { + + CurrentElem currelem(iel,myproc,Level,VV,eqn,*mesh,eqn->GetMLProb().GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,eqn->GetMLProb().GetQrule(currelem.GetDim())); @@ -124,19 +136,8 @@ double ComputeIntegral (const uint Level, const MultiLevelMeshTwo* mesh, const S xyz_refbox._FEord = MESH_ORDER; xyz_refbox._ndof = mymsh->el->GetElementDofNumber(ZERO_ELEM,BIQUADR_FE); xyz_refbox.Allocate(); - - - double integral = 0.; - - const uint el_ngauss = eqn->GetMLProb().GetQrule(currelem.GetDim()).GetGaussPointsNumber(); - -//parallel sum - const uint nel_e = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level+1]; - const uint nel_b = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level]; - for (uint iel=0; iel < (nel_e - nel_b); iel++) { - - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -154,7 +155,7 @@ double ComputeIntegral (const uint Level, const MultiLevelMeshTwo* mesh, const S if ( Tdes._eqnptr != NULL ) Tdes.GetElemDofs(); else Tdes._qtyptr->FunctionDof(Tdes,0.,&xyz_refbox._val_dofs[0]); - + const uint el_ngauss = eqn->GetMLProb().GetQrule(currelem.GetDim()).GetGaussPointsNumber(); for (uint qp = 0; qp < el_ngauss; qp++) { @@ -228,7 +229,18 @@ double ComputeNormControl (const uint Level, const MultiLevelMeshTwo* mesh, cons Mesh *mymsh = eqn->GetMLProb()._ml_msh->GetLevel(Level); - CurrentElem currelem(Level,VV,eqn,*mesh,eqn->GetMLProb().GetElemType()); + + + double integral = 0.; + + +//parallel sum + const uint nel_e = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level+1]; + const uint nel_b = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level]; + + for (int iel=0; iel < (nel_e - nel_b); iel++) { + + CurrentElem currelem(iel,myproc,Level,VV,eqn,*mesh,eqn->GetMLProb().GetElemType()); currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,eqn->GetMLProb().GetQrule(currelem.GetDim())); @@ -254,19 +266,10 @@ double ComputeNormControl (const uint Level, const MultiLevelMeshTwo* mesh, cons xyz_refbox._ndof = mymsh->el->GetElementDofNumber(ZERO_ELEM,BIQUADR_FE); xyz_refbox.Allocate(); - - double integral = 0.; - //loop over the geom el types - const uint el_ngauss = eqn->GetMLProb().GetQrule(currelem.GetDim()).GetGaussPointsNumber(); - -//parallel sum - const uint nel_e = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level+1]; - const uint nel_b = mesh->_off_el[mesh_vb][mesh->_NoLevels*myproc+Level]; - - for (int iel=0; iel < (nel_e - nel_b); iel++) { - - currelem.SetDofobjConnCoords(myproc,iel); + const uint el_ngauss = eqn->GetMLProb().GetQrule(currelem.GetDim()).GetGaussPointsNumber(); + + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); diff --git a/applications/OptimalControl/tempopt/TempQuantities.cpp b/applications/OptimalControl/tempopt/TempQuantities.cpp index 7b88164c5..6a4fd381f 100644 --- a/applications/OptimalControl/tempopt/TempQuantities.cpp +++ b/applications/OptimalControl/tempopt/TempQuantities.cpp @@ -54,7 +54,6 @@ Velocity::Velocity(std::string name_in, QuantityMap& qtymap_in, uint dim_in, uin //============================================================= -///analytical velocity for Hartmann flow // difference between get_par and optsys: // in both cases you are "dynamic" somehow @@ -112,34 +111,6 @@ void Velocity::Function_txyz(const double /*t*/,const double* xp, double* func) } - -//============================================================= -void Velocity::strain_txyz(const double /*t*/, const double* xyz,double strain[][DIMENSION]) const { - -//here, tau is a tensor, so tau dot n is a vector which in general has a NORMAL and a TANGENTIAL component - - const double Lref = _qtymap.GetInputParser()->get("Lref"); - double ILref = 1./Lref; - const double lye = _qtymap.GetMeshTwo()->GetDomain()->_domain_rtmap.get("lye"); -// const double x=xyz[0]; - const double y=xyz[1]; - if (_qtymap.GetMeshTwo()->get_dim() == 3) { - const double z=xyz[2]; - } - - strain[0][0] = 0.; //ux,x - strain[0][1] = strain[1][0] = 0. ;//0.5*(uy,x+ux,y) - strain[1][1] = 0.*(-(lye*ILref-y)); //uy,y - if (_qtymap.GetMeshTwo()->get_dim() == 3) { - strain[0][2] = strain[2][0] = 0. ; //0.5*(uz,x+ux,z) - strain[1][2] = strain[2][1] = 0. ; //0.5*(uy,z+uz,y) - strain[2][2] = 0. ; //uz,z - } - -return; -} - - //============================================================= /// prescribed pressure at the boundary //no initial condition for pressure is required, because it has no time derivative @@ -259,27 +230,6 @@ void Temperature::Function_txyz(const double/* t*/, const double* xp,double* tem } -// ================================================= - //the coordinates (x,y,z,t) of the VOLUME domain are NON-dimensional - //and the function value must be nondimensional as well - //-----Nonhomogeneous Neumann------- - // Qflux = - k grad(T) by definition -// QfluxDOTn>0: energy flows outside (cooling) QfluxDOTn<0: energy flows inside (heating) -void Temperature::heatflux_txyz(const double /*t*/, const double* /*xyz*/, double* qflux) const { - -// std::cout << "Temperature: Heatflux, check which coordinates are passed in here" << std::endl; -// Box* box= static_cast(_qtymap._phys._mesh->GetDomain()); -// const double thetaz = box->_domain_rtmap.get("thetaz"); - - qflux[0]=-2.1*0./**cos(thetaz)*/; - qflux[1]=0./**sin(thetaz)*/; - if (_qtymap.GetMeshTwo()->get_dim() == 3) { - qflux[2]=0.; - } - - return; - } - void Velocity::bc_flag_txyz(const double t, const double* xp, std::vector & bc_flag) const { diff --git a/applications/OptimalControl/tempopt/TempQuantities.hpp b/applications/OptimalControl/tempopt/TempQuantities.hpp index 66ce30698..9135ecafc 100644 --- a/applications/OptimalControl/tempopt/TempQuantities.hpp +++ b/applications/OptimalControl/tempopt/TempQuantities.hpp @@ -22,11 +22,6 @@ class Temperature : public Quantity { void bc_flag_txyz(const double t, const double* xp, std::vector & flag) const; void initialize_xyz(const double* xp, std::vector & value) const; -//specific function - //this is the function of the IMPOSED DERIVATIVE of TEMPERATURE, aka heat flux - void heatflux_txyz(const double t,const double* xyz, double* qflux) const; - - }; //=============== @@ -98,13 +93,6 @@ class Velocity : public Quantity { void bc_flag_txyz(const double t, const double* xp, std::vector & flag) const; void initialize_xyz(const double* xp, std::vector & value) const; - //specific function - //this is the STRAIN DERIVATIVE of VELOCITY, so it must stay here - //from the physical and also mathematical point of view - //the shape funcs of the same order as v will then be used and so on -//multi-dimensional arrays must have bounds for all dimensions except the first - void strain_txyz(const double t, const double* xp,double strain[][DIMENSION]) const; //TODO convert this double array - }; diff --git a/applications/OptimalControl/tempopt/main.cpp b/applications/OptimalControl/tempopt/main.cpp index a5a70ccc4..ec366715a 100644 --- a/applications/OptimalControl/tempopt/main.cpp +++ b/applications/OptimalControl/tempopt/main.cpp @@ -135,7 +135,17 @@ void GenMatRhsNS(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gr ml_msh.SetDomain(&mybox); MultiLevelSolution ml_sol(&ml_msh); - ml_sol.AddSolution("FAKE",LAGRANGE,SECOND,0); + ml_sol.AddSolution("Qty_Temperature",LAGRANGE,SECOND,0); + ml_sol.AddSolution("Qty_TempLift",LAGRANGE,SECOND,0); + ml_sol.AddSolution("Qty_TempAdj",LAGRANGE,SECOND,0); + ml_sol.AddSolutionVector(ml_msh.GetDimension(),"Qty_Velocity",LAGRANGE,SECOND,0); + ml_sol.AddSolution("Qty_Pressure",LAGRANGE,FIRST,0); + ml_sol.AddSolution("Qty_TempDes",LAGRANGE,SECOND,0,false); //this is not going to be an Unknown! + + ml_sol.Initialize("All"); /// @todo you have to call this before you can print @todo I can also call it after instantiation MLProblem + ml_sol.SetWriter(VTK); + std::vector print_vars(1); print_vars[0] = "All"; // we should find a way to make this easier + ml_sol.GetWriter()->write(files.GetOutputPath(),"biquadratic",print_vars); MultiLevelProblem ml_prob(&ml_sol); ml_prob.SetMeshTwo(&mesh); @@ -150,13 +160,13 @@ void GenMatRhsNS(MultiLevelProblem &ml_prob, unsigned Level, const unsigned &gr // not all the Quantities need to be unknowns of an equation SystemTwo & eqnNS = ml_prob.add_system("Eqn_NS"); - eqnNS.AddSolutionToSystemPDE("FAKE"); +// eqnNS.AddSolutionToSystemPDE("FAKE"); //do the vector version of it eqnNS.AddUnknownToSystemPDE(&velocity); eqnNS.AddUnknownToSystemPDE(&pressure); eqnNS.SetAssembleFunction(GenMatRhsNS); SystemTwo & eqnT = ml_prob.add_system("Eqn_T"); - eqnT.AddSolutionToSystemPDE("FAKE"); +// eqnT.AddSolutionToSystemPDE("FAKE"); eqnT.AddUnknownToSystemPDE(&temperature); eqnT.AddUnknownToSystemPDE(&templift); eqnT.AddUnknownToSystemPDE(&tempadj);//the order in which you add defines the order in the matrix as well, so it is in tune with the assemble function diff --git a/doc/devel_rules.txt b/doc/devel_rules.txt index 83eed58d6..e19b64d73 100644 --- a/doc/devel_rules.txt +++ b/doc/devel_rules.txt @@ -96,4 +96,36 @@ git config --global user.email "name.surname@example.com" The maintainers will decide what to do with the pull request and possibly it will be merged to master Periodically, sync the master in the fork with the master in the main femus repository +################################## + +WORKFLOW for updating the master in the FORK from the master in the MAIN REPO + +In github: + +Go in the MAIN REPO + +Click on "Pull requests" (on the right) +Click on "New pull request" (green button) +Click on "compare across forks" +The "base fork" is going to be the FORK master branch +The "head fork" is going to be the MAIN master branch +Click on "create pull request" +Add some title for it +Make sure that the branches can be AUTOMATICALLY MERGED (otherwise you have to solve the conflicts using command line...) +Click on "Merge pull request" (you find it by scrolling towards the bottom of the page) +Click on "Confirm merge" + +From command line: + +To be added + +################################## + +How to install HDFVIEW + +Follow the instruction in the HDFVIEW website + +At the end of the build process, you should have an hdfview.sh script inside the bin directory of your build. +You must open that script and change the INSTALLDIR variable with the path of this build (by default, /usr/local is in it) +Another alternative to view the content of an HDF5 file is to use the 'h5dump' utility shipped with HDF5 (installed through PETSc, for instance) diff --git a/src/algebra/LinearEquation.cpp b/src/algebra/LinearEquation.cpp index afaface63..fa402fc6c 100644 --- a/src/algebra/LinearEquation.cpp +++ b/src/algebra/LinearEquation.cpp @@ -42,8 +42,6 @@ LinearEquation::LinearEquation(Mesh *other_msh){ _EPSC = NULL; _RES = NULL; _RESC = NULL; - _PP = NULL; - _RR = NULL; _KK = NULL; _CC = NULL; } @@ -237,14 +235,6 @@ void LinearEquation::DeletePde() { delete _CC; } - - if (_msh->GetLevel()>0) { - if(_PP) - delete _PP; - if(_RR) - delete _RR; - } - if(_EPS) delete _EPS; diff --git a/src/algebra/LinearEquation.hpp b/src/algebra/LinearEquation.hpp index 50bc24a49..8576faad3 100644 --- a/src/algebra/LinearEquation.hpp +++ b/src/algebra/LinearEquation.hpp @@ -79,7 +79,7 @@ class LinearEquation : public ParallelObject { // member data Mesh *_msh; NumericVector *_EPS, *_EPSC, *_RES, *_RESC; - SparseMatrix *_KK, *_PP,*_RR, *_CC; + SparseMatrix *_KK, *_CC; vector < vector > KKoffset; vector < unsigned > KKghostsize; vector < vector < int> > KKghost_nd; diff --git a/src/enums/FETypeEnum.hpp b/src/enums/FETypeEnum.hpp index e95386c7e..76971ef97 100644 --- a/src/enums/FETypeEnum.hpp +++ b/src/enums/FETypeEnum.hpp @@ -8,9 +8,7 @@ enum FEType {QQ=0, LL, KK }; #define QL 3 #define BIQUADR_FE 2 -#define BIQUADR_TYPEEL 3 #define LINEAR_FE 0 -#define LINEAR_TYPEEL 0 #define NFE_FAMS 5 diff --git a/src/equations/BoundaryConditions.cpp b/src/equations/BoundaryConditions.cpp index 1d120fb15..6c83d13b1 100644 --- a/src/equations/BoundaryConditions.cpp +++ b/src/equations/BoundaryConditions.cpp @@ -223,16 +223,18 @@ void BoundaryConditions::GenerateBdc() { for (uint Level=0; Level <_dofmap->_mesh._NoLevels;Level++) { //loop over the levels Mesh *mymsh = _dofmap->_eqn->GetMLProb()._ml_msh->GetLevel(Level); - CurrentElem currelem(Level,BB,_dofmap->_eqn,_dofmap->_mesh,_dofmap->_eqn->GetMLProb().GetElemType()); - currelem.SetMesh(mymsh); const uint el_nnodes_b = mymsh->el->GetElementFaceDofNumber(ZERO_ELEM,ZERO_FACE,BIQUADR_FE); for (uint isubd=0; isubd<_dofmap->_mesh._NoSubdom; ++isubd) { uint iel_b = _dofmap->_mesh._off_el[BB][ _dofmap->_mesh._NoLevels*isubd + Level]; uint iel_e = _dofmap->_mesh._off_el[BB][ _dofmap->_mesh._NoLevels*isubd + Level+1]; + for (uint iel=0; iel < (iel_e - iel_b); iel++) { + + CurrentElem currelem(iel,isubd,Level,BB,_dofmap->_eqn,_dofmap->_mesh,_dofmap->_eqn->GetMLProb().GetElemType()); + currelem.SetMesh(mymsh); - currelem.SetDofobjConnCoords(isubd,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); for (uint ivar=0; ivar< _dofmap->_n_vars; ivar++) bc_flag[ivar] = DEFAULT_BC_FLAG; //this is necessary here to re-clean! diff --git a/src/equations/CurrentElem.cpp b/src/equations/CurrentElem.cpp index 0ce609e57..d634aeb2a 100644 --- a/src/equations/CurrentElem.cpp +++ b/src/equations/CurrentElem.cpp @@ -27,13 +27,15 @@ namespace femus { - CurrentElem::CurrentElem(const uint level, const uint vb, const SystemTwo * eqn_in, const MultiLevelMeshTwo& mesh, const std::vector< std::vector > & elem_type_in ): + CurrentElem::CurrentElem(const uint iel_in, const uint iproc_in, const uint level, const uint vb, const SystemTwo * eqn_in, const MultiLevelMeshTwo& mesh, const std::vector< std::vector > & elem_type_in ): _eqn(eqn_in), _mesh(mesh), _dim(_mesh.get_dim()-vb), _elem_type(elem_type_in[mesh.get_dim()-vb -1]), _mesh_vb(vb), - _Level(level) + _Level(level), + _iel(iel_in), + _proc(iproc_in) { //========== Current "Geometric Element" ======================== @@ -227,14 +229,14 @@ void CurrentElem::PrintOrientation() const { } // ===================================================================================== - void CurrentElem::SetDofobjConnCoords(const uint isubd_in,const uint iel) { + void CurrentElem::SetDofobjConnCoords() { const uint mydim = _mesh.get_dim(); const uint el_nnodes = _el_conn.size(); for (uint n=0; nGetNumberOfElements(); + +// elem * myel = _mesh_new->el; +// unsigned kel = _mesh_new->IS_Mts2Gmt_elem[_iel]; +// short unsigned kelt = myel->GetElementType(kel); + + +// unsigned order_ind = ml_sol->GetSolutionType(ml_sol->GetIndex("U")); +// unsigned nve = myel->GetElementDofNumber(kel,order_ind); +// +// +// +// _el_conn_new.resize(nve); +// +// for (unsigned i=0;iGetElementVertexIndex(kel,i)-1u):(kel+i*nel); +// unsigned inode=myel->GetElementVertexIndex(kel,i)-1u; + +// // dof metis +// /*metis_node2*/_el_conn_new[i] = mymsh->GetMetisDof(inode,BIQUADR_FE); +// metis_node1[i] = mymsh->GetMetisDof(inode,SolType[2*dim]); + + +// dofsVAR[j+dim][i]= mylsyspde->GetKKDof(indVAR[j+dim],indexVAR[j+dim],inode); + +// } +// +// +// +// return; +// } + + + // ===================================================================================== +// void CurrentElem::SetCoords_two() { +// +// for (unsigned i=0;iGetElementVertexIndex(kel,i)-1u; +// dof metis +// unsigned inode_Metis=mymsh->GetMetisDof(inode,2); +// metis_node2[i]=inode_Metis; +// +// unsigned inode_Metis=mymsh->GetMetisDof(inode,2); +// flag to know if the node "inode" lays on the fluid-solid interface +// solidmark[i]=myel->GetNodeRegion(inode); // to check +// for(int j=0; j_coordinate->_Sol[j])(inode_Metis) + (*mysolution->_Sol[indVAR[j]])(inode_Metis); +// Old coordinates (Moving frame) +// vx_old[j][i]= (*mymsh->_coordinate->_Sol[j])(inode_Metis) + (*mysolution->_SolOld[indVAR[j]])(inode_Metis); +// Fixed coordinates (Reference frame) +// vx_hat[j][i]= (*mymsh->_coordinate->_Sol[j])(inode_Metis); +// displacement dofs +// dofsVAR[j][i]= mylsyspde->GetKKDof(indVAR[j],indexVAR[j],inode); +// velocity dofs +// dofsVAR[j+dim][i]= mylsyspde->GetKKDof(indVAR[j+dim],indexVAR[j+dim],inode); +// } +// } +// +// +// +// +// return; +// } + + + } //end namespace femus diff --git a/src/equations/CurrentElem.hpp b/src/equations/CurrentElem.hpp index 5da41c98b..457e03e25 100644 --- a/src/equations/CurrentElem.hpp +++ b/src/equations/CurrentElem.hpp @@ -36,7 +36,7 @@ class CurrentQuantity; public: - CurrentElem(const uint level, const uint vb, const SystemTwo*, const MultiLevelMeshTwo& mesh, const std::vector< std::vector > & elem_type); + CurrentElem(const uint iel_in, const uint iproc_in, const uint level, const uint vb, const SystemTwo*, const MultiLevelMeshTwo& mesh, const std::vector< std::vector > & elem_type); inline const uint GetVb() const { return _mesh.get_dim() - _dim; @@ -81,7 +81,7 @@ class CurrentQuantity; return _KeM; } - void SetDofobjConnCoords(const uint isubd_in,const uint iel); + void SetDofobjConnCoords(); void SetMidpoint(); @@ -117,23 +117,26 @@ class CurrentQuantity; // ======================================================================================== //========== Current "EQUATION" Element (ql are TOGETHER ): needs the EQUATION ======================== + uint _el_n_dofs; DenseMatrix _KeM; DenseVector _FeM; - uint _el_n_dofs; std::vector _el_dof_indices; // this must become a vect of vect std::vector _bc_eldofs; // this must become a vect of vect // ======================================================================================== //========== Current Geometric Element: needs the MESH ======================== std::vector _el_conn; /// vector of the global nodes for that element [NNDS]; - std::vector _el_conn_new; uint _vol_iel_DofObj; /// i need to put the element also. + std::vector _el_conn_new; std::vector _xx_nds; /// vector of the node coordinates for that element [_spacedimension*NNDS]; // this must become a vect of vect std::vector _el_xm; /// element center point [_spacedimension]; const uint _dim; //spatial dimension of the current element (can be different from the mesh dimension!) const uint _mesh_vb; //index for the mesh const uint _Level; //the level to which the element belongs + + const uint _iel; //the index of the element (input for the parallel partition) + const uint _proc; }; @@ -172,11 +175,6 @@ class CurrentQuantity; //in fact this is used for many element loop, but just to retrieve the geometrical properties //like coords, middle point, etc. -//Ok, so far we have VB, so if we want only BB or VV we should actually make a vector -// of two functions like this -//TODO i should TEMPLATIZE over VB!!! YES IT IS TRUE, IN FACT BASICALLY EVERYTHING -//DEPENDS ON VB, and when you use one you dont use the other!!! - //======================= //ok here we would need a "REFERENCE REAL ELEMENT" //the current element contains the absolute coordinates, @@ -188,6 +186,9 @@ class CurrentQuantity; //The curr geometric is basically filled with the MESH class //The curr fe is basically filled with the EQUATION class - +//==================================================== +// The differences with the other applications are that: +// the element matrix is split into blocks, while here it is only one +// Therefore, also the el_dof_indices and the bc_eldof flags are split as vector of vectors #endif \ No newline at end of file diff --git a/src/equations/CurrentQuantity.cpp b/src/equations/CurrentQuantity.cpp index 384fc5e9b..13e91448c 100644 --- a/src/equations/CurrentQuantity.cpp +++ b/src/equations/CurrentQuantity.cpp @@ -290,6 +290,25 @@ void CurrentQuantity::GetElemDofs() { } +// void CurrentQuantity::GetElemDofs_two() { +// +// unsigned SolType2 = ml_sol->GetSolutionType(ml_sol->GetIndex(_qtyptr->_name)); +// unsigned nve2 = _currEl._mesh_new->el->GetElementDofNumber(kel,SolType2); +// +// +// for (unsigned i=0;iel->GetMeshDof(kel,i,SolType2); +// unsigned inode_Metis=_currEl._mesh_new->GetMetisDof(inode,SolType2); +// for(int j=0; j<_currEl._mesh_new->GetDimension(); j++) { +// // velocity dofs +// Soli[indexVAR[j]][i] = (*mysolution->_Sol[indVAR[j]])(inode_Metis); +// // dofsVAR[j][i] = myLinEqSolver->GetKKDof(indVAR[j],indexVAR[j],inode); +// } +// } + + + + // clearly _el_average must be allocated already!!! void CurrentQuantity::SetElemAverage() { diff --git a/src/equations/LinearImplicitSystem.cpp b/src/equations/LinearImplicitSystem.cpp index b3b5e4351..55bf29a6d 100644 --- a/src/equations/LinearImplicitSystem.cpp +++ b/src/equations/LinearImplicitSystem.cpp @@ -109,6 +109,8 @@ void LinearImplicitSystem::init() { AddVariableToBeSolved("All"); } + + /// @deprecated // this function is like init but it doesn't call InitPDE void LinearImplicitSystem::init_two() { @@ -124,6 +126,15 @@ void LinearImplicitSystem::init_two() { // _LinSolver[i]->InitPde(_SolSystemPdeIndex,_ml_sol->GetSolType(), // _ml_sol->GetSolName(),&_solution[i]->_Bdc,_gridr,_gridn,_SparsityPattern); // } + + + _PP.resize(_gridn); + _RR.resize(_gridn); + for(unsigned i=0;i<_gridn;i++){ + _PP[i]=NULL; + _RR[i]=NULL; + } + // // for (unsigned ig=1; ig<_gridn; ig++) { // BuildProlongatorMatrix(ig); diff --git a/src/equations/LinearImplicitSystem.hpp b/src/equations/LinearImplicitSystem.hpp index af5ea8fd5..9afea01f2 100644 --- a/src/equations/LinearImplicitSystem.hpp +++ b/src/equations/LinearImplicitSystem.hpp @@ -168,9 +168,10 @@ class LinearImplicitSystem : public ImplicitSystem { /** enforce sparcity pattern for setting uncoupled variables and save on memory allocation **/ void SetSparsityPattern(vector < bool > other_sparcity_pattern); + vector < SparseMatrix* > _PP, _RR; /// @todo put it back to protected + protected: - vector < SparseMatrix* > _PP, _RR; /** Create the Prolongator matrix for the Multigrid solver */ void Prolongator(const unsigned &gridf); diff --git a/src/equations/SystemTwo.cpp b/src/equations/SystemTwo.cpp index 2e452fe47..c9e267faa 100644 --- a/src/equations/SystemTwo.cpp +++ b/src/equations/SystemTwo.cpp @@ -229,17 +229,19 @@ void SystemTwo::Initialize() { for (uint Level = 0; Level< GetGridn(); Level++) { - Mesh *mymsh = GetMLProb()._ml_msh->GetLevel(Level); - CurrentElem currelem(Level,VV,this,GetMLProb().GetMeshTwo(),GetMLProb().GetElemType()); - currelem.SetMesh(mymsh); - const uint el_dof_objs = NVE[ GetMLProb().GetMeshTwo()._geomelem_flag[currelem.GetDim()-1] ][BIQUADR_FE]; + Mesh *mymsh = GetMLProb()._ml_msh->GetLevel(Level); + const unsigned myproc = mymsh->processor_id(); uint iel_b = GetMLProb().GetMeshTwo()._off_el[VV][ GetMLProb().GetMeshTwo()._iproc*GetGridn() + Level ]; uint iel_e = GetMLProb().GetMeshTwo()._off_el[VV][ GetMLProb().GetMeshTwo()._iproc*GetGridn() + Level + 1]; for (uint iel=0; iel < (iel_e - iel_b); iel++) { - currelem.SetDofobjConnCoords(GetMLProb().GetMeshTwo()._iproc,iel); + CurrentElem currelem(iel,myproc,Level,VV,this,GetMLProb().GetMeshTwo(),GetMLProb().GetElemType()); + currelem.SetMesh(mymsh); + const uint el_dof_objs = NVE[ GetMLProb().GetMeshTwo()._geomelem_flag[currelem.GetDim()-1] ][BIQUADR_FE]; + + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); for (uint q=0; q < _UnknownQuantitiesVector.size() ; q++) { diff --git a/src/mesh/GambitIO.cpp b/src/mesh/GambitIO.cpp index a0a297887..45b16094f 100644 --- a/src/mesh/GambitIO.cpp +++ b/src/mesh/GambitIO.cpp @@ -24,7 +24,7 @@ namespace femus { - const unsigned GambitIO::GambitToFemusVertexIndex[6][27]= + const unsigned GambitIO::GambitToFemusVertexIndex[N_GEOM_ELS][27]= { { 4,16,0,15,23,11,7,19,3, @@ -46,7 +46,7 @@ namespace femus { }; -const unsigned GambitIO::GambitToFemusFaceIndex[6][6]= +const unsigned GambitIO::GambitToFemusFaceIndex[N_GEOM_ELS][6]= { {0,4,2,5,3,1}, {0,1,2,3}, diff --git a/src/mesh/GambitIO.hpp b/src/mesh/GambitIO.hpp index 4e5608133..28ec75567 100644 --- a/src/mesh/GambitIO.hpp +++ b/src/mesh/GambitIO.hpp @@ -19,6 +19,7 @@ // Local includes #include "MeshInput.hpp" +#include "GeomElTypeEnum.hpp" namespace femus { @@ -53,10 +54,10 @@ class GambitIO : public MeshInput private: /** Map from Gambit vertex index to Femus vertex index */ - static const unsigned GambitToFemusVertexIndex[6][27]; + static const unsigned GambitToFemusVertexIndex[N_GEOM_ELS][27]; /** Map from Gambit face index to Femus face index */ - static const unsigned GambitToFemusFaceIndex[6][6]; + static const unsigned GambitToFemusFaceIndex[N_GEOM_ELS][6]; }; diff --git a/src/mesh/SalomeIO.cpp b/src/mesh/SalomeIO.cpp index 09ad633c1..759bd69c0 100644 --- a/src/mesh/SalomeIO.cpp +++ b/src/mesh/SalomeIO.cpp @@ -34,30 +34,66 @@ namespace femus { const uint SalomeIO::max_length = 100; ///@todo this length of the menu string is conservative enough... + //How to determine a general connectivity: + //you have to align the element with respect to the x-y-z (or xi-eta-zeta) reference frame, + //and then look at the order in the med file. + //For every node there is a location, and you have to put that index in that x-y-z location. + //Look NOT at the NUMBERING, but at the ORDER! - const unsigned SalomeIO::SalomeToFemusVertexIndex[6][27]= +// SALOME HEX27 +// 1------17-------5 +// /| /| +// / | / | +// 8 | 21 12 | +// / 9 22 / 13 +// / | / | +// 0------16-------4 | +// | 20 | 26 | 25 | +// | 2------18-|-----6 zeta +// | / | / ^ +// 11 / 24 15 / | eta +// | 10 23 | 14 | / +// | / | / | / +// |/ |/ |/ +// 3-------19------7 -------> xi + +// TEMPLATE HEX27 (for future uses) +// X------X--------X +// /| /| +// / | / | +// X | X X | +// / X X / X +// / | / | +// X------X--------X | +// | X | X | X | +// | X------X--|-----X zeta +// | / | / ^ +// X / X X / | eta +// | X X | X | / +// | / | / | / +// |/ |/ |/ +// X-------X-------X -------> xi + + + const unsigned SalomeIO::SalomeToFemusVertexIndex[N_GEOM_ELS][27]= { - { - 4,16,0,15,23,11,7,19,3, - 12,20,8,25,26,24,14,22,10, - 5,17,1,13,21,9,6,18,2 - }, - { - 0,4,1,6,5, - 2,7,8,9,3 - }, + + {4,7,3,0,5,6,2,1,15,19,11,16,13,18,9,17,12,14,10,8,23,25,22,24,20,21,26}, //HEX27 + {0,1,2,3,4,5,6,7,8,9}, //TET10 + { 3, 11,5, 9, 10,4, 12,17,14,15,16,13, 0, 8, 2, 6, 7, 1 - }, - {0,4,1,5,2,6,3,7,8}, - {0,3,1,4,2,5}, - {0,2,1} + }, //WEDGE18 + + {0,1,2,3,4,5,6,7,8}, //QUAD9 + {0,1,2,3,4,5}, //TRI6 + {0,1,2} //EDGE3 }; -const unsigned SalomeIO::SalomeToFemusFaceIndex[6][6]= +const unsigned SalomeIO::SalomeToFemusFaceIndex[N_GEOM_ELS][6]= { {0,4,2,5,3,1}, {0,1,2,3}, @@ -193,9 +229,9 @@ void SalomeIO::read(const std::string& name, vector < vector < double> > &coords // SET NUMBER OF ELEMENTS mesh.SetNumberOfElements(n_elements); - int *conn_map5 = new int[dim_conn]; + int *conn_map = new int[dim_conn]; std::cout << " Number of elements in med file " << n_elements << std::endl; - status=H5Dread(dtset,H5T_NATIVE_INT,H5S_ALL,H5S_ALL,H5P_DEFAULT,conn_map5); + status=H5Dread(dtset,H5T_NATIVE_INT,H5S_ALL,H5S_ALL,H5P_DEFAULT,conn_map); if(status !=0) {std::cout << "SalomeIO::read: connectivity not found"; abort();} H5Dclose(dtset); @@ -209,28 +245,28 @@ void SalomeIO::read(const std::string& name, vector < vector < double> > &coords if (nve==27) { type_elem_flag[0]=type_elem_flag[3]=true; mesh.el->AddToElementNumber(1,"Hex"); - mesh.el->SetElementType(iel,0); + mesh.el->SetElementType(iel,HEX); } else if (nve==10) { type_elem_flag[1]=type_elem_flag[4]=true; mesh.el->AddToElementNumber(1,"Tet"); - mesh.el->SetElementType(iel,1); + mesh.el->SetElementType(iel,TET); } else if (nve==18) { type_elem_flag[2]=type_elem_flag[3]=type_elem_flag[4]=true; mesh.el->AddToElementNumber(1,"Wedge"); - mesh.el->SetElementType(iel,2); + mesh.el->SetElementType(iel,WEDGE); } else if (nve==9) { type_elem_flag[3]=true; mesh.el->AddToElementNumber(1,"Quad"); - mesh.el->SetElementType(iel,3); + mesh.el->SetElementType(iel,QUAD); } else if (nve==6 && mesh.GetDimension()==2) { type_elem_flag[4]=true; mesh.el->AddToElementNumber(1,"Triangle"); - mesh.el->SetElementType(iel,4); + mesh.el->SetElementType(iel,TRI); } else if (nve==3 && mesh.GetDimension()==1) { mesh.el->AddToElementNumber(1,"Line"); - mesh.el->SetElementType(iel,5); + mesh.el->SetElementType(iel,LINE); } else { std::cout<<"Error! Invalid element type in reading File!"< > &coords } for (unsigned i=0; iGetElementType(iel)][i]; - mesh.el->SetElementVertexIndex(iel,inode,conn_map5[iel+i*n_elements]); + mesh.el->SetElementVertexIndex(iel,inode,conn_map[iel+i*n_elements]); } } - -// // // // // // // Adding Connectivity to mesh structure (Libmesh) -// // // // // mesh.reserve_elem(n_elements); -// // // // // // read the elements -// // // // // for(int iel=0; iel<(int)n_elements; ++iel) { -// // // // // // add the elements to the mesh -// // // // // Elem* elem = mesh.add_elem(Elem::build(eletype.type).release()); -// // // // // // add node pointers to the elements -// // // // // for(int i=0; iset_node(eletype.nodes[i])= mesh.node_ptr(conn_map[Node_el*iel+i]); -// // // // // elem->set_node(eletype.nodes[i])= mesh.node_ptr(conn_map5[iel+i*n_elements]-1); -// // // // // } -// // // // // } // // end read ELEMENT/CELL **************** B @@ -271,7 +294,7 @@ void SalomeIO::read(const std::string& name, vector < vector < double> > &coords // clean - delete [] conn_map5; + delete [] conn_map; } @@ -345,7 +368,7 @@ void SalomeIO::read(const std::string& name, vector < vector < double> > &coords // unsigned iel,iface; // inf>>iel>>str2>>iface; // iel--; -// iface=SalomeIO::GambitToFemusFaceIndex[mesh.el->GetElementType(iel)][iface-1u]; +// iface=SalomeIO::SalomeToFemusFaceIndex[mesh.el->GetElementType(iel)][iface-1u]; // mesh.el->SetFaceElementIndex(iel,iface,value); // } // inf >> str2; diff --git a/src/mesh/SalomeIO.hpp b/src/mesh/SalomeIO.hpp index 97410af89..740c611b8 100644 --- a/src/mesh/SalomeIO.hpp +++ b/src/mesh/SalomeIO.hpp @@ -19,6 +19,8 @@ // Local includes #include "MeshInput.hpp" +#include "GeomElTypeEnum.hpp" + #ifdef HAVE_HDF5 #include "hdf5.h" #endif @@ -56,10 +58,10 @@ class SalomeIO : public MeshInput private: /** Map from Salome vertex index to Femus vertex index */ - static const unsigned SalomeToFemusVertexIndex[6][27]; + static const unsigned SalomeToFemusVertexIndex[N_GEOM_ELS][27]; /** Map from Salome face index to Femus face index */ - static const unsigned SalomeToFemusFaceIndex[6][6]; + static const unsigned SalomeToFemusFaceIndex[N_GEOM_ELS][6]; /** Determine mesh dimension from mesh file */ void FindDimension(hid_t gid, const std::string menu_name,hsize_t n_fem_type); diff --git a/src/meshGencase/GenCase.cpp b/src/meshGencase/GenCase.cpp index 60d437b17..9de00bc4e 100644 --- a/src/meshGencase/GenCase.cpp +++ b/src/meshGencase/GenCase.cpp @@ -2115,7 +2115,7 @@ void GenCase::ComputeMaxElXNode() { // ========================================= -void GenCase::ReadMGOps(const std::string output_path, const SystemTwo * mysys) { +void GenCase::ReadMGOps(const std::string output_path, SystemTwo * mysys) { std::string f_matrix = DEFAULT_F_MATRIX; std::string f_rest = DEFAULT_F_REST; @@ -2297,7 +2297,7 @@ void GenCase::ReadMGOps(const std::string output_path, const SystemTwo * mysys) // So, if the node_dof was not filled correctly, then when you -void GenCase::ReadMatrix(const std::string& namefile, const SystemTwo * mysys) { +void GenCase::ReadMatrix(const std::string& namefile, SystemTwo * mysys) { for (uint Level = 0; Level< mysys->GetGridn(); Level++) { @@ -2490,7 +2490,7 @@ void GenCase::ReadMatrix(const std::string& namefile, const SystemTwo * mysys) //============================= //This function depends on _iproc -void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) { +void GenCase::ReadProl(const std::string& name, SystemTwo * mysys) { for (uint Level = 1; Level< mysys->GetGridn(); Level++) { @@ -2573,7 +2573,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) { uint off_proc = mysys->GetMLProb().GetMeshTwo()._iproc*mysys->GetGridn(); - mysys->_LinSolver[Lev_f]->_PP = SparseMatrix::build().release(); + mysys->_PP[Lev_f] = SparseMatrix::build().release(); // // // _Prl[ Lev_f ]->init(0,0,0,0); //TODO BACK TO A REASONABLE INIT // local matrix dimension @@ -2650,7 +2650,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) { std::cout << "Printing Prolongator ===========" << std::endl; pattern.print(); - mysys->_LinSolver[Lev_f]->_PP->update_sparsity_pattern_old(pattern); + mysys->_PP[Lev_f]->update_sparsity_pattern_old(pattern); //=========== VALUES =================== DenseMatrix *valmat; @@ -2670,7 +2670,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) { for (uint j=0; j_LinSolver[Lev_f]->_PP->add_matrix(*valmat,tmp,ind); + mysys->_PP[Lev_f]->add_matrix(*valmat,tmp,ind); delete valmat; } } @@ -2691,7 +2691,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) { pattern.clear(); - mysys->_LinSolver[Lev_f]->_PP->close(); + mysys->_PP[Lev_f]->close(); // if (mysys->GetMLProb().GetMeshTwo()._iproc==0) _Prl[ Lev_f ]->print_personal(); // _Prl[ Lev_f ]->print_graphic(false); //TODO should pass this true or false as a parameter } //end levels @@ -2727,7 +2727,7 @@ void GenCase::ReadProl(const std::string& name, const SystemTwo * mysys) { //AAA fai molta attenzione: per esplorare la node_dof devi usare Lev_c e Lev_f, //perche' sono legati ai DOF (devi pensare che la questione del mesh e' gia' risolta) -void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) { +void GenCase::ReadRest(const std::string& name, SystemTwo * mysys) { for (uint Level = 0; Level< mysys->GetGridn() - 1; Level++) { @@ -2809,7 +2809,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) { uint off_proc=mysys->GetGridn()*mysys->GetMLProb().GetMeshTwo()._iproc; - mysys->_LinSolver[Lev_c]->_RR = SparseMatrix::build().release(); + mysys->_RR[Lev_c] = SparseMatrix::build().release(); // // // _Rst[Lev_c]->init(0,0,0,0); //TODO BACK TO A REASONABLE INIT //we have to do this before appropriately!!! int nrowt=0;int nclnt=0; @@ -2873,7 +2873,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) { std::cout << "Printing Restrictor ===========" << std::endl; pattern.print(); - mysys->_LinSolver[Lev_c]->_RR->update_sparsity_pattern_old(pattern); //TODO see + mysys->_RR[Lev_c]->update_sparsity_pattern_old(pattern); //TODO see // _Rst[Lev_c]->close(); // if (mysys->GetMLProb().GetMeshTwo()._iproc==0) _Rst[Lev_c]->print_personal(); //there is no print function for rectangular matrices, and print_personal doesnt seem to be working... // la print stampa il contenuto, ma io voglio solo stampare lo sparsity pattern! @@ -2901,7 +2901,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) { for (uint i1=0;i1_bcond._bc[irow_top]*Rest_val[fe][ j+len[fe][i] ]; - mysys->_LinSolver[Lev_c]->_RR->add_matrix(*valmat,tmp,ind); + mysys->_RR[Lev_c]->add_matrix(*valmat,tmp,ind); delete valmat; }// end dof loop } // end var loop @@ -2921,7 +2921,7 @@ void GenCase::ReadRest(const std::string& name, const SystemTwo * mysys) { pattern.clear(); - mysys->_LinSolver[Lev_c]->_RR->close(); + mysys->_RR[Lev_c]->close(); // if (mysys->GetMLProb().GetMeshTwo()._iproc==0) _Rst[Lev_c]->print_personal(std::cout); // _Rst[Lev_c]->print_graphic(false); // TODO should pass this true or false as a parameter diff --git a/src/meshGencase/GenCase.hpp b/src/meshGencase/GenCase.hpp index d7ba36da4..157bf9f21 100644 --- a/src/meshGencase/GenCase.hpp +++ b/src/meshGencase/GenCase.hpp @@ -47,10 +47,10 @@ class GenCase : public MultiLevelMeshTwo { void ComputeAndPrintProl(const std::string output_path); void ComputeAndPrintRest(const std::string output_path); - static void ReadMGOps(const std::string output_path, const SystemTwo * mysys); - static void ReadMatrix(const std::string& name, const SystemTwo * mysys); - static void ReadProl(const std::string& name, const SystemTwo * mysys); - static void ReadRest(const std::string& name, const SystemTwo * mysys); + static void ReadMGOps(const std::string output_path, SystemTwo * mysys); + static void ReadMatrix(const std::string& name, SystemTwo * mysys); + static void ReadProl(const std::string& name, SystemTwo * mysys); + static void ReadRest(const std::string& name, SystemTwo * mysys); void CreateMeshStructuresLevSubd(const std::string output_path); void Delete(); diff --git a/src/solution/MultiLevelSolution.cpp b/src/solution/MultiLevelSolution.cpp index e5e52b5b9..f2609fe7a 100644 --- a/src/solution/MultiLevelSolution.cpp +++ b/src/solution/MultiLevelSolution.cpp @@ -131,6 +131,18 @@ void MultiLevelSolution::AddSolution(const char name[], const FEFamily fefamily, } } + void MultiLevelSolution::AddSolutionVector(const unsigned n_components, const std::string name, const FEFamily fefamily, const FEOrder order, unsigned tmorder, const bool &Pde_type) { + + for (unsigned i=0; i _solution; + /** Flag to tell whether the BC function has been set */ + bool _bdc_func_set; + /** This group of vectors has the size of the number of added solutions */ vector< vector > _boundaryconditions; vector< vector > _ishomogeneous; vector< vector > _nonhomogeneousbcfunction; - bool _bdc_func_set; - unsigned short _gridn; - vector _SolType; + vector _SolType; /* Tells the FE index */ vector _family; vector _order; vector _SolName; vector _BdcType; vector _SolTmorder; - vector _PdeType; + vector _PdeType; /*Tells whether the Solution is an unknown of a PDE or not*/ vector _TestIfPressure; vector _SolPairIndex; diff --git a/src/solution/VTKWriter.cpp b/src/solution/VTKWriter.cpp index 6dc49290b..6ed8067e0 100644 --- a/src/solution/VTKWriter.cpp +++ b/src/solution/VTKWriter.cpp @@ -125,7 +125,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const short unsigned ielt = _ml_mesh->GetLevel(ig)->el->GetElementType(kel); for (unsigned j=0; j<_ml_mesh->GetLevel(ig)->el->GetElementDofNumber(kel,index); j++) { counter++; - unsigned loc_vtk_conn = map_pr[j]; + unsigned loc_vtk_conn = FemusToVTKorToXDMFConn[j]; unsigned jnode=_ml_mesh->GetLevel(ig)->el->GetMeshDof(kel, loc_vtk_conn, index); unsigned jnodeMetis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode, index); if( jnodeMetis < offset_iprc ){ // check if jnodeMetis is a ghost node @@ -220,7 +220,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const if ( ig == _gridn-1u || 0 == _ml_mesh->GetLevel(ig)->el->GetRefinedElementIndex(kel)) { short unsigned ielt=_ml_mesh->GetLevel(ig)->el->GetElementType(kel); for (unsigned j=0; j<_ml_mesh->GetLevel(ig)->el->GetElementDofNumber(kel,index); j++) { - unsigned loc_vtk_conn = map_pr[j]; + unsigned loc_vtk_conn = FemusToVTKorToXDMFConn[j]; unsigned jnode=_ml_mesh->GetLevel(ig)->el->GetMeshDof(kel, loc_vtk_conn, index); unsigned jnodeMetis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode, index); if( jnodeMetis < offset_iprc ){ @@ -241,7 +241,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const if ( ig == _gridn-1u || 0 == _ml_mesh->GetLevel(ig)->el->GetRefinedElementIndex(kel)) { short unsigned ielt=_ml_mesh->GetLevel(ig)->el->GetElementType(kel); for (unsigned j=0; j<_ml_mesh->GetLevel(ig)->el->GetElementDofNumber(kel,index); j++) { - unsigned loc_vtk_conn = map_pr[j]; + unsigned loc_vtk_conn = FemusToVTKorToXDMFConn[j]; unsigned jnode=_ml_mesh->GetLevel(ig)->el->GetMeshDof(kel, loc_vtk_conn, index); unsigned jnodeMetis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode, index); if( jnodeMetis < offset_iprc ){ @@ -294,7 +294,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const unsigned kel = _ml_mesh->GetLevel(ig)->IS_Mts2Gmt_elem[iel]; if ( ig == _gridn-1u || 0 == _ml_mesh->GetLevel(ig)->el->GetRefinedElementIndex(kel)) { for (unsigned j=0; j<_ml_mesh->GetLevel(ig)->el->GetElementDofNumber(kel,index); j++) { - unsigned loc_vtk_conn = map_pr[j]; + unsigned loc_vtk_conn = FemusToVTKorToXDMFConn[j]; unsigned jnode=_ml_mesh->GetLevel(ig)->el->GetMeshDof(kel, loc_vtk_conn, index); unsigned jnodeMetis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode, index); var_conn[icount] = (jnodeMetis >= offset_iprc )? jnodeMetis - offset_iprc + offset_nvt : nvt0 + (ghost_counter++); @@ -566,7 +566,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const if ( ig == _gridn-1u || 0 == _ml_mesh->GetLevel(ig)->el->GetRefinedElementIndex(kel)) { short unsigned ielt=_ml_mesh->GetLevel(ig)->el->GetElementType(kel); for (unsigned j=0; j<_ml_mesh->GetLevel(ig)->el->GetElementDofNumber(kel,index); j++) { - unsigned loc_vtk_conn = map_pr[j]; + unsigned loc_vtk_conn = FemusToVTKorToXDMFConn[j]; unsigned jnode=_ml_mesh->GetLevel(ig)->el->GetMeshDof(kel, loc_vtk_conn, index); unsigned jnodeMetis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode, index); if( jnodeMetis < offset_iprc ){ @@ -802,7 +802,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const // for (unsigned iel=0; iel<_ml_mesh->GetLevel(ig)->GetNumberOfElements(); iel++) { // if (ig==_gridn-1u || _ml_mesh->GetLevel(ig)->el->GetRefinedElementIndex(iel)==0) { // for (unsigned j=0; j<_ml_mesh->GetLevel(ig)->el->GetElementDofNumber(iel,index); j++) { -// unsigned loc_vtk_conn = map_pr[j]; +// unsigned loc_vtk_conn = FemusToVTKorToXDMFConn[j]; // unsigned jnode=_ml_mesh->GetLevel(ig)->el->GetElementVertexIndex(iel,loc_vtk_conn)-1u; // unsigned jnodeMetis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode,index); // var_conn[icount] = offset_nvt + jnodeMetis; diff --git a/src/solution/Writer.cpp b/src/solution/Writer.cpp index bccccfd99..1adb6ae2f 100644 --- a/src/solution/Writer.cpp +++ b/src/solution/Writer.cpp @@ -2,7 +2,7 @@ Program: FEMUS Module: Writer - Authors: Eugenio Aulisa, Simone Bnà + Authors: Eugenio Aulisa, Simone Bnà, Giorgio Bornia Copyright (c) FEMTTU All rights reserved. @@ -33,7 +33,9 @@ namespace femus { - + + const unsigned Writer::FemusToVTKorToXDMFConn[27] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,21,20,22,24,25,26}; + Writer::Writer( MultiLevelSolution* ml_sol ): _ml_sol(ml_sol), _ml_mesh(ml_sol->_ml_msh) { diff --git a/src/solution/Writer.hpp b/src/solution/Writer.hpp index 6e55b6961..52e0eb774 100644 --- a/src/solution/Writer.hpp +++ b/src/solution/Writer.hpp @@ -2,7 +2,7 @@ Program: FEMUS Module: Writer - Authors: Eugenio Aulisa, Simone Bnà + Authors: Eugenio Aulisa, Simone Bnà, Giorgio Bornia Copyright (c) FEMTTU All rights reserved. @@ -27,9 +27,6 @@ namespace femus { - // map from our connectivity to vtk-connectivity for paraview visualization //TODO move this to the appropriate place - const unsigned map_pr[27] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,21,20,22,24,25,26}; - //------------------------------------------------------------------------------ // Forward declarations //------------------------------------------------------------------------------ @@ -82,6 +79,8 @@ namespace femus { int _gridr; + /** map from femus connectivity to vtk-connectivity for paraview visualization */ + static const unsigned FemusToVTKorToXDMFConn[27]; private: diff --git a/src/solution/XDMFWriter.cpp b/src/solution/XDMFWriter.cpp index 10271d2e2..1288d68a0 100644 --- a/src/solution/XDMFWriter.cpp +++ b/src/solution/XDMFWriter.cpp @@ -42,10 +42,11 @@ namespace femus { - const std::string XDMFWriter::type_el[4][6] = {{"Hexahedron","Tetrahedron","Wedge","Quadrilateral","Triangle","Polyline"}, - {"Hexahedron_20","Tetrahedron_10","Not_implemented","Quadrilateral_8","Triangle_6","Edge_3"}, - {"Not_implemented","Not_implemented","Not_implemented","Not_implemented","Not_implemented","Not_implemented"}, - {"Hexahedron_27","Not_implemented","Not_implemented","Quadrilateral_9","Triangle_6","Edge_3"}}; + const std::string XDMFWriter::type_el[3][N_GEOM_ELS] = {{"Hexahedron","Tetrahedron","Wedge","Quadrilateral","Triangle","Polyline"}, //linear + {"Hexahedron_20","Tetrahedron_10","Not_implemented","Quadrilateral_8","Triangle_6","Edge_3"}, //serendipity + {"Hexahedron_27","Tetrahedron_10","Not_implemented","Quadrilateral_9","Triangle_6","Edge_3"}}; //tensor-product quadratic (some real, some fake) + /// @todo Tri6 and Tet10 are actually SERENDIPITY, not TENSOR-PRODUCT QUADRATIC. + // The corresponding tensor-product ones should be Tri7 and Tet14 const std::string XDMFWriter::_nodes_name = "/NODES"; const std::string XDMFWriter::_elems_name = "/ELEMS"; @@ -68,25 +69,21 @@ void XDMFWriter::write(const std::string output_path, const char order[], const print_all += !(vars[ivar].compare("All")) + !(vars[ivar].compare("all")) + !(vars[ivar].compare("ALL")); } - unsigned index=0; unsigned index_nd=0; if(!strcmp(order,"linear")) { //linear - index=0; index_nd=0; } else if(!strcmp(order,"quadratic")) { //quadratic - index=1; index_nd=1; } - else if(!strcmp(order,"biquadratic")) { //biquadratic - index=3; + else if(!strcmp(order,"biquadratic")) { //tensor-product quadratic (real and fake) index_nd=2; } /// @todo I assume that the mesh is not mixed std::string type_elem; unsigned elemtype = _ml_mesh->GetLevel(_gridn-1u)->el->GetElementType(ZERO_ELEM); - type_elem = XDMFWriter::type_el[index][elemtype]; + type_elem = XDMFWriter::type_el[index_nd][elemtype]; if (type_elem.compare("Not_implemented") == 0) { @@ -108,7 +105,7 @@ void XDMFWriter::write(const std::string output_path, const char order[], const nel+=_ml_mesh->GetLevel(_gridn-1u)->GetNumberOfElements(); unsigned icount; - unsigned el_dof_number = _ml_mesh->GetLevel(_gridn-1u)->el->GetElementDofNumber(0,index_nd); + unsigned el_dof_number = _ml_mesh->GetLevel(_gridn-1u)->el->GetElementDofNumber(ZERO_ELEM,index_nd); int * var_conn = new int [nel*el_dof_number]; std::vector< int > var_proc(nel); float *var_el_f = new float [nel]; @@ -277,7 +274,7 @@ void XDMFWriter::write(const std::string output_path, const char order[], const if ( ig == _gridn-1u || _ml_mesh->GetLevel(ig)->el->GetRefinedElementIndex(iel) == 0) { int ndofs = _ml_mesh->GetLevel(ig)->el->GetElementDofNumber(iel,index_nd); for (unsigned j = 0; j < ndofs; j++) { - unsigned vtk_loc_conn = map_pr[j]; + unsigned vtk_loc_conn = FemusToVTKorToXDMFConn[j]; unsigned jnode = _ml_mesh->GetLevel(ig)->el->GetElementVertexIndex(iel,vtk_loc_conn)-1u; unsigned jnode_Metis = _ml_mesh->GetLevel(ig)->GetMetisDof(jnode,index_nd); var_conn[icount] = offset_conn + jnode_Metis; @@ -1487,9 +1484,9 @@ void XDMFWriter::PrintXDMFTopGeom(std::ofstream& out, #ifdef HAVE_HDF5 - uint n_children, order_typeel; - if (order_fe == BIQUADR_FE) { n_children = 1; order_typeel = BIQUADR_TYPEEL;} - else if (order_fe == LINEAR_FE) {n_children = NRE[mesh._eltype_flag[vb]]; order_typeel = LINEAR_TYPEEL;} + uint n_children; + if (order_fe == BIQUADR_FE) { n_children = 1; } + else if (order_fe == LINEAR_FE) {n_children = NRE[mesh._eltype_flag[vb]]; } else { std::cout << "Mesh Not supported" << std::endl; abort(); } uint nel = mesh._n_elements_vb_lev[vb][Level]; @@ -1498,7 +1495,7 @@ void XDMFWriter::PrintXDMFTopGeom(std::ofstream& out, std::ostringstream coord_lev; coord_lev << "_L" << Level; - PrintXDMFTopology(out,top_file.str(),hdf_field.str(), type_el[order_typeel][mesh._eltype_flag[vb]], nel*n_children, nel*n_children, NVE[mesh._eltype_flag[vb]][order_fe]); + PrintXDMFTopology(out,top_file.str(),hdf_field.str(), type_el[order_fe][mesh._eltype_flag[vb]], nel*n_children, nel*n_children, NVE[mesh._eltype_flag[vb]][order_fe]); PrintXDMFGeometry(out,geom_file.str(),_nodes_name+"/COORD/X",coord_lev.str(),"X_Y_Z","Float",mesh._NoNodesXLev[Level],1); diff --git a/src/solution/XDMFWriter.hpp b/src/solution/XDMFWriter.hpp index 819a018cb..ef61adb9c 100644 --- a/src/solution/XDMFWriter.hpp +++ b/src/solution/XDMFWriter.hpp @@ -146,7 +146,7 @@ class XDMFWriter : public Writer { private: - static const std::string type_el[4][6]; + static const std::string type_el[3][N_GEOM_ELS]; static const std::string _nodes_name; static const std::string _elems_name; diff --git a/src/utils/Math.hpp b/src/utils/Math.hpp index 43d7f1f88..92ba94e25 100644 --- a/src/utils/Math.hpp +++ b/src/utils/Math.hpp @@ -107,10 +107,22 @@ const uint myproc = ml_prob.GetMeshTwo()._iproc; Mesh *mymsh = ml_prob._ml_msh->GetLevel(Level); - CurrentElem currelem(Level,vb,NULL,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); //element without equation + + double integral = 0.; + + +//parallel sum + const uint nel_e = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level+1]; + const uint nel_b = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level]; + + for (uint iel=0; iel < (nel_e - nel_b); iel++) { + + CurrentElem currelem(iel,myproc,Level,vb,NULL,ml_prob.GetMeshTwo(),ml_prob.GetElemType()); //element without equation currelem.SetMesh(mymsh); CurrentGaussPointBase & currgp = CurrentGaussPointBase::build(currelem,ml_prob.GetQrule(currelem.GetDim())); + const uint el_ngauss = ml_prob.GetQrule(currelem.GetDim()).GetGaussPointsNumber(); + //========= DOMAIN MAPPING CurrentQuantity xyz(currgp); xyz._dim = ml_prob.GetMeshTwo().get_dim(); @@ -118,18 +130,7 @@ const uint myproc = ml_prob.GetMeshTwo()._iproc; xyz._ndof = currelem.GetElemType(xyz._FEord)->GetNDofs(); xyz.Allocate(); - double integral = 0.; - -//loop over the geom el types - const uint el_ngauss = ml_prob.GetQrule(currelem.GetDim()).GetGaussPointsNumber(); - -//parallel sum - const uint nel_e = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level+1]; - const uint nel_b = ml_prob.GetMeshTwo()._off_el[mesh_vb][ml_prob.GetMeshTwo()._NoLevels*myproc+Level]; - - for (uint iel=0; iel < (nel_e - nel_b); iel++) { - - currelem.SetDofobjConnCoords(myproc,iel); + currelem.SetDofobjConnCoords(); currelem.SetMidpoint(); currelem.ConvertElemCoordsToMappingOrd(xyz); @@ -160,9 +161,9 @@ double myval_g = pt2func(time,xyz._val_g); std::cout << std::endl << " ^^^^^^^^^^^^^^^^^L'integrale sul processore "<< myproc << " vale: " << integral << std::endl; - double weights_sum = 0.; - for (uint qp = 0; qp < el_ngauss; qp++) weights_sum += ml_prob.GetQrule(currelem.GetDim()).GetGaussWeight(qp); - std::cout << std::endl << " ^^^^^^^^^^^^^^^^^ La somma dei pesi vale: " << weights_sum << std::endl; +// double weights_sum = 0.; +// for (uint qp = 0; qp < el_ngauss; qp++) weights_sum += ml_prob.GetQrule(currelem.GetDim()).GetGaussWeight(qp); +// std::cout << std::endl << " ^^^^^^^^^^^^^^^^^ La somma dei pesi vale: " << weights_sum << std::endl; double J=0.; #ifdef HAVE_MPI diff --git a/unittests/testSalomeIO/input/OneEdge3.med b/unittests/testSalomeIO/input/OneEdge3.med new file mode 100644 index 000000000..8d88efcdb Binary files /dev/null and b/unittests/testSalomeIO/input/OneEdge3.med differ diff --git a/unittests/testSalomeIO/input/OneHex27.med b/unittests/testSalomeIO/input/OneHex27.med index 1dc3ca571..b63a6b695 100644 Binary files a/unittests/testSalomeIO/input/OneHex27.med and b/unittests/testSalomeIO/input/OneHex27.med differ diff --git a/unittests/testSalomeIO/input/OneTet10.med b/unittests/testSalomeIO/input/OneTet10.med new file mode 100644 index 000000000..9962fe441 Binary files /dev/null and b/unittests/testSalomeIO/input/OneTet10.med differ diff --git a/unittests/testSalomeIO/input/OneTri6.med b/unittests/testSalomeIO/input/OneTri6.med new file mode 100644 index 000000000..b5071648c Binary files /dev/null and b/unittests/testSalomeIO/input/OneTri6.med differ diff --git a/unittests/testSalomeIO/input/StudyHex27.hdf b/unittests/testSalomeIO/input/StudyHex27.hdf index 8e804e29b..b5fd1111c 100644 Binary files a/unittests/testSalomeIO/input/StudyHex27.hdf and b/unittests/testSalomeIO/input/StudyHex27.hdf differ diff --git a/unittests/testSalomeIO/input/StudyTet10.hdf b/unittests/testSalomeIO/input/StudyTet10.hdf new file mode 100644 index 000000000..8e04f3735 Binary files /dev/null and b/unittests/testSalomeIO/input/StudyTet10.hdf differ diff --git a/unittests/testSalomeIO/input/StudyTri6.hdf b/unittests/testSalomeIO/input/StudyTri6.hdf new file mode 100644 index 000000000..7283f5293 Binary files /dev/null and b/unittests/testSalomeIO/input/StudyTri6.hdf differ diff --git a/unittests/testSalomeIO/testSalomeIO.cpp b/unittests/testSalomeIO/testSalomeIO.cpp index 7f348bff8..f143257b6 100644 --- a/unittests/testSalomeIO/testSalomeIO.cpp +++ b/unittests/testSalomeIO/testSalomeIO.cpp @@ -21,7 +21,7 @@ int main(int argc,char **args) { double Lref = 1.; MultiLevelMesh ml_msh; - ml_msh.ReadCoarseMesh(infile.c_str(),"seventh",Lref); + ml_msh.ReadCoarseMesh(infile.c_str(),"fifth",Lref); ml_msh.SetWriter(XDMF); ml_msh.GetWriter()->write(DEFAULT_OUTPUTDIR,"biquadratic");