Merge pull request #5 from FeMTTU/master

Merge from main repo hex27

applications/FSI/include/FSIassembly.hpp

@@ -163,7 +163,7 @@ void AssembleMatrixResFSI(MultiLevelProblem &ml_prob, unsigned level, const unsi
   vector <unsigned> indexVAR(2*dim+1);
   //vector <unsigned> indCOORD(dim);
   vector <unsigned> indVAR(3*dim+1);  
-  vector <unsigned> SolType(3*dim+1);  
+  vector <unsigned> SolType(3*dim+1);  /// @todo is this really 3*d + 1 or 2*d+1 would be enough?
 
   for(unsigned ivar=0; ivar<dim; ivar++) {
     //indCOORD[ivar]=ml_prob.GetIndex(&coordname[ivar][0]);

applications/OptimalControl/fe_test/EqnT.cpp

@@ -65,7 +65,7 @@
 
   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()));
 
@@ -106,7 +106,7 @@
     currelem.Mat().zero();
     currelem.Rhs().zero(); 
 
-    currelem.SetDofobjConnCoords(myproc,iel);
+    currelem.SetDofobjConnCoords();
     currelem.SetMidpoint();
 
     currelem.ConvertElemCoordsToMappingOrd(xyz);

applications/OptimalControl/fe_test/TempQuantities.cpp

@@ -41,26 +41,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<Box*>(_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<int> & bc_flag) const  {
 // T' and its adjoint must be Dirichlet homogeneous everywhere on the boundary, by definition.

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<int> & flag) const;
   void initialize_xyz(const double* xp, std::vector<double> & 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;
-
 
 };
 

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);

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);

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);

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);

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);

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);

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<Box*>(_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;
-  }
 
 
 

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<int> & flag) const;
   void initialize_xyz(const double* xp, std::vector<double> & 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;
-
-
 };