add pbc minimal constraint ele

NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp

@@ -360,11 +360,12 @@ nonLinearGeneralPeriodicBC::nonLinearGeneralPeriodicBC(const nonLinearGeneralPer
     
 void nonLinearGeneralPeriodicBC::setPBCNormal(const double n1, const double n2, const double n3){
     _pbcNormal[0] = n1;
-    _pbcNormal[1] = n1;
-    _pbcNormal[2] = n2;
+    _pbcNormal[1] = n2;
+    _pbcNormal[2] = n3;
     if (_pbcNormal.norm() >0){
       _pbcNormal.normalize();
     }
+    printf("set PBC normal = [%f %f %f]\n",_pbcNormal[0],_pbcNormal[1],_pbcNormal[2]);
 };
 
 nonLinearDisplacementBC::nonLinearDisplacementBC( const int tag, const int dim):

NonLinearSolver/periodicBC/directionalConstraintElement.cpp

@@ -514,6 +514,7 @@ void directionalPBCSupplementConstraint::getLinearConstraints(std::map<Dof,DofAf
   getLinearConstraints(g, con);
 };
 
+
 directionalPBCLagrangeConstraintElement::directionalPBCLagrangeConstraintElement(nonLinearMicroBC* mbc, FunctionSpaceBase* space, FunctionSpaceBase* mspace, 
 															MVertex* v1, std::vector<MVertex*>& vlist, const SVector3& dir):
 															constraintElement(mbc,-1), _periodicSpace(space),_multSpace(mspace),

NonLinearSolver/periodicBC/pbcConstraintElement.cpp

@@ -469,6 +469,9 @@ lagrangeConstraintElement::lagrangeConstraintElement(nonLinearMicroBC* mbc, Func
 	}
 	
 	getPBCKinematicMatrix(_L,_XX,_S);
+  
+  //this->print();
+  //Msg::Info("project point %d is %f %f %f",_vp->getNum(),point[0],point[1],point[2]);
 };
 
 

NonLinearSolver/periodicBC/pbcCreateConstraints.cpp

@@ -1728,7 +1728,62 @@ void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
 	else{
 		Msg::Fatal("createShiftedLagrangeConstraintElementGroup has not been implemented for 3D problems");
 	}
+};
+
+void pbcConstraintElementGroup::createGeneralPBCLagrangeConstraintElementGroup(){
+  nonLinearGeneralPeriodicBC* gpbc = static_cast<nonLinearGeneralPeriodicBC*>(_solver->getMicroBC());
+  const SVector3& pbcNormal = gpbc->getPBCNormal();
+  int dim = _solver->getMicroBC()->getDim();
+  if (dim == 2){
+    SPoint3 p1 = _v1->point();
+		SPoint3 p2 = _v2->point();
+		SPoint3 p3 = _v3->point();
+    
+    SPoint3 p = planeDirProject(p3,p1,p2,pbcNormal);
+    Msg::Info("found p1 = [%f %f %f], p2 = [%f %f %f], p = [%f %f %f]",p1[0],p1[1],p1[2],p2[0],p2[1],p2[2],p[0],p[1],p[2]);
+    
+    SVector3 vecp2p(p2,p);
+    SVector3 vecp2p1(p2,p1);
+    if (dot(vecp2p,vecp2p1) < 0){
+      // base is constraint following l12, l23
+      printf("interpolation based following l12,l23\n");
+      std::set<MVertex*> newVirVertices;
+      InterpolationOperations::createPolynomialBasePoints(gpbc,_v2,_v1,_v3,_v5,l12,l23,l15,_xBase,_yBase,_zBase,newVirVertices);
+      _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
+      
+      // form condition
+      this->lagrangeFormCondition(l12.begin(),l12.end(),_xBase);
+      this->lagrangeFormCondition(l23.begin(),l23.end(),_yBase);
+      
+      // periodic BC
+      std::set<MVertex*> allPositiveVertex;
+      allPositiveVertex.insert(l41.begin(),l41.end());
+      allPositiveVertex.insert(l34.begin(),l34.end());
+      this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,pbcNormal,true);
+      
+    }
+    else{
+      // base is constraint following l12, l41
+      printf("interpolation based following l12,l41\n");
+      std::set<MVertex*> newVirVertices;
+      InterpolationOperations::createPolynomialBasePoints(gpbc,_v1,_v2,_v4,_v5,l12,l41,l15,_xBase,_yBase,_zBase,newVirVertices);
+      _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
+      
+      // form condition
+      this->lagrangeFormCondition(l12.begin(),l12.end(),_xBase);
+      this->lagrangeFormCondition(l41.begin(),l41.end(),_yBase);
+      
+      // periodic BC
+      std::set<MVertex*> allPositiveVertex;
+      allPositiveVertex.insert(l23.begin(),l23.end());
+      allPositiveVertex.insert(l34.begin(),l34.end());
+      this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,pbcNormal,true);
+    }
   }
+  else {
+    Msg::Fatal("createShiftedLagrangeConstraintElementGroup has not been implemented for 3D problems");
+  }
+};
 
 void pbcConstraintElementGroup::createLagrangeConstraintElementGroup(){
   int dim = _solver->getMicroBC()->getDim();
@@ -1796,10 +1851,28 @@ void pbcConstraintElementGroup::createLinearDisplacementConstraintElementGroup()
 void pbcConstraintElementGroup::createMinimalKinematicConstraintElementGroup(){
   const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
   int size = bgroup.size();
+  /*
   for (int i=0; i<size; i++){
 		constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,-1,bgroup[i]);
 		_allConstraint.insert(ppc);
   }
+   */
+  for (int i=0; i< size/2; i++){
+    groupOfElements* gPlus = bgroup[i+size/2];
+    groupOfElements* gMinus = bgroup[i];
+
+    if (i == 0){
+      constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,-1,gPlus);
+      _allConstraint.insert(ppc);
+      ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,-1,gMinus);
+      _allConstraint.insert(ppc);
+    }
+    else{
+       constraintElement* ppc = new periodicSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,-1,gPlus,gMinus);
+       _allConstraint.insert(ppc);
+       
+    }
+  }
 };
 
 void pbcConstraintElementGroup::createLinearDisplacementConstraintElementGroup_DG(){
@@ -2529,6 +2602,53 @@ void pbcConstraintElementGroup::createPBCConstraintGroup(){
   this->createCornerConstraintElementGroup();
 };
 
+void pbcConstraintElementGroup::createGeneralPBCConstraintGroup(){
+  nonLinearGeneralPeriodicBC* gpbc = static_cast<nonLinearGeneralPeriodicBC*>(_solver->getMicroBC());
+  const SVector3& pbcNormal = gpbc->getPBCNormal();
+  bool considerUsualPBC = false;
+  if (pbcNormal.norm() < 1e-10){
+    printf("because pbc vector is null, usual PBC is used\n");
+    considerUsualPBC = true;
+  }
+  
+  if (!considerUsualPBC){
+    SPoint3 p1 = _v1->point();
+		SPoint3 p2 = _v2->point();
+		SPoint3 p3 = _v3->point();
+    SPoint3 p4 = _v4->point();
+    
+    SVector3 p1p2(p1,p2);
+    SVector3 p1p4(p1,p4);
+    
+    if ( (fabs(dot(p1p2,pbcNormal))< 1e-6) or (fabs(dot(p1p4,pbcNormal))< 1e-6)){
+      considerUsualPBC = true;
+      printf("because pbc vector coincide with boundary normal, usual PBC is used\n");
+    }
+  }
+  
+  
+  if (considerUsualPBC){
+   this->createPBCConstraintGroup();
+  }
+  else{
+    printf("general PBC is used with Lagrange interpolation method\n");
+    gpbc->setPBCMethod(nonLinearPeriodicBC::LIM);
+    
+    if (gpbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
+      Msg::Info("Imposing PBC by lagrange formulation");
+      this->createGeneralPBCLagrangeConstraintElementGroup();
+    }
+    else{
+      Msg::Fatal("this method is not implemented to impose shiftedPBC");
+    };
+    
+    // add minimal kinematic
+    printf("add minimal kinematic BC \n");
+    createMinimalKinematicConstraintElementGroup();
+  }
+  
+};
+
 void pbcConstraintElementGroup::createShiftedPBCConstraintGroup(){
   nonLinearShiftedPeriodicBC* spbc = static_cast<nonLinearShiftedPeriodicBC*>(_solver->getMicroBC());
   const SVector3& normal = spbc->getShiftPBCNormal();
@@ -2580,6 +2700,9 @@ void pbcConstraintElementGroup::createConstraintGroup(){
 		}
     else if (_solver->getMicroBC()->getType() == nonLinearMicroBC::ShiftedPBC){
       this->createShiftedPBCConstraintGroup();
+    }
+    else if (_solver->getMicroBC()->getType() == nonLinearMicroBC::GeneralPBC){
+      this->createGeneralPBCConstraintGroup();
     }
 		else if (_solver->getMicroBC()->getType()==nonLinearMicroBC::LDBC){
 			Msg::Info("Imposing linear displacement BC");

NonLinearSolver/periodicBC/pbcCreateConstraints.h

@@ -240,6 +240,7 @@ class pbcConstraintElementGroup{
     void createCubicSplineConstraintElementGroup();
     void createLagrangeConstraintElementGroup();
 		void createShiftedLagrangeConstraintElementGroup();
+    void createGeneralPBCLagrangeConstraintElementGroup();
     void createFEConstraintElementGroup();
     void createLinearDisplacementConstraintElementGroup();
 		void createProjectConstraintElementGroup();
@@ -258,6 +259,7 @@ class pbcConstraintElementGroup{
     
     void createPBCConstraintGroup();
     void createShiftedPBCConstraintGroup();
+    void createGeneralPBCConstraintGroup();
 
     SVector3 getUniformDisp(MVertex* v);
     SVector3 getPerturbation(dofManager<double>* pmanager, MVertex* v);

NonLinearSolver/periodicBC/pbcSupplementConstraint.cpp

@@ -234,6 +234,196 @@ void supplementConstraint::getLinearConstraints(std::map<Dof,DofAffineConstraint
   getLinearConstraints(g, con);
 };
 
+
+periodicSupplementConstraint::periodicSupplementConstraint(nonLinearMicroBC* mbc, FunctionSpaceBase* space, FunctionSpaceBase* mspace,
+												const int c, groupOfElements* gplus, groupOfElements* gminus):
+									supplementConstraintBase(mbc, space,mspace,c),_gPlus(gplus),_gMinus(gminus){
+  std::vector<MVertex*> vPlus, vMinus;
+  fullVector<double> weightMinus,weightPlus;
+                    
+	supplementConstraint::computeWeight(_gPlus,vPlus,weightPlus);
+  supplementConstraint::computeWeight(_gMinus,vMinus,weightMinus);
+  
+  int sizeMinus = vMinus.size();
+  int sizePlus = vPlus.size();
+  int sizever = sizePlus+sizeMinus;
+  
+  _v.resize(sizever);
+  _weight.resize(sizever);
+  for (int i=0; i< sizePlus; i++){
+    _v[i] = vPlus[i];
+    _weight(i) = weightPlus(i);
+  }
+  for (int i=0; i< sizeMinus; i++){
+    _v[i+sizePlus] = vMinus[i];
+    _weight(i+sizePlus) = -weightMinus(i);
+  }
+	
+  // to obatain a good positive vertex
+  std::vector<Dof> oneDofPerVertex;
+  for (int i=0; i< _v.size(); i++){
+    std::vector<Dof> keys;
+    getKeysFromVertex(_periodicSpace,_v[i],getComp(),keys);
+    oneDofPerVertex.push_back(keys[0]);
+  }
+  _positive = -1;
+	double maxVal = 0.;
+	bool found = false;
+  for (int i=0; i<sizever; i++){
+    if ((fabs(_weight(i))>maxVal) && (constraintElement::allPositiveDof.find(oneDofPerVertex[i]) == constraintElement::allPositiveDof.end())){
+      _positive = i;
+			maxVal = fabs(_weight(i));
+			found = true;
+    }
+  };
+	if (!found){
+		Msg::Fatal("All Dof are numerated as positive dof in other constraints supplementConstraint::supplementConstraint");
+	}
+
+  // add positive to positive vertex list
+  std::vector<Dof> Keys;
+  getKeysFromVertex(_periodicSpace,_v[_positive],getComp(),Keys);
+  for (int ik=0; ik < Keys.size(); ik++){
+    if (constraintElement::allPositiveDof.find(Keys[ik]) == constraintElement::allPositiveDof.end()){
+      constraintElement::allPositiveDof.insert(Keys[ik]);
+    }
+    else{
+      Msg::Warning("Dof on vertex was chosen as positive one in other constraint element:");
+    }
+  }
+
+  double invs = 1./(_weight(_positive));
+  _weight.scale(invs);
+  _Xmean*= 0;
+  for (int i=0; i<_v.size();i++){
+    MVertex* v = _v[i];
+    SPoint3 pnt = v->point();
+    for( int j=0; j<3; j++){
+      _Xmean[j] += _weight(i)*pnt[j];
+    }
+  }
+	
+	if (_mbc->getOrder() == 2){
+		_XXmean*= 0;
+		for (int i=0; i<_v.size();i++){
+			MVertex* v = _v[i];
+			SPoint3 pnt = v->point();
+			for( int j=0; j<3; j++){
+				for (int k=0; k<3; k++){
+					_XXmean(j,k) += _weight(i)*0.5*pnt[j]*pnt[k];
+				}
+			}
+		}
+	}
+	
+  // add a virtual vertex for identify the constraint
+  _tempVertex = new MVertex(_v[_positive]->x(),_v[_positive]->y(),_v[_positive]->z());
+				
+	_C.resize(_ndofs,sizever*_ndofs);
+	_C.setAll(0.);
+
+	_Cbc.resize(_ndofs, (sizever-1)*_ndofs);
+	_Cbc.setAll(0.);
+
+	int col = 0;
+	int colbc = 0;
+	for (int i=0; i<_weight.size(); i++){
+		for (int j=0; j<_ndofs; j++){
+			_C(j,col) = _weight(i);
+			if (i != _positive){
+				_Cbc(j,colbc) = -1.*_weight(i);
+				colbc++;
+			}
+			col++;
+		}
+	}
+	
+	getPointKinematicMatrix(_Xmean,_XXmean,_S);
+};
+
+
+void periodicSupplementConstraint::getConstraintKeys(std::vector<Dof>& key) const {
+  for (int i=0; i<_v.size(); i++)
+      getKeysFromVertex(_periodicSpace,_v[i],getComp(),key);
+}
+void periodicSupplementConstraint::getMultiplierKeys(std::vector<Dof>& key) const{
+  getKeysFromVertex(_multSpace,_tempVertex,getComp(),key);
+}
+void periodicSupplementConstraint::getConstraintMatrix(fullMatrix<double>& m) const {
+  m = _C;
+};		// matrix C
+void periodicSupplementConstraint::getDependentMatrix(fullMatrix<double>& m) const{
+  m = _Cbc;
+}
+
+// for constraint form  u = C u* + Sb*FM + M uc
+void periodicSupplementConstraint::getDependentKeys(std::vector<Dof>& keys) const {
+  getKeysFromVertex(_periodicSpace,_v[_positive],getComp(),keys);
+}; // left real dofs
+
+void periodicSupplementConstraint::getIndependentKeys(std::vector<Dof>& keys) const{
+  for (int i=0; i<_v.size(); i++)
+    if (i!= _positive)
+      getKeysFromVertex(_periodicSpace,_v[i],getComp(),keys);
+}
+
+
+
+void periodicSupplementConstraint::getLinearConstraints(const fullVector<double>& g,
+                   std::map<Dof,DofAffineConstraint<double> >& con) const
+{
+
+  std::vector<std::vector<Dof> > k;
+  for (int i=0; i< _v.size(); i++){
+    std::vector<Dof> ktemp;
+    getKeysFromVertex(_periodicSpace,_v[i],getComp(),ktemp);
+    k.push_back(ktemp);
+  };
+
+  std::vector<Dof> kp;
+  this->getDependentKeys(kp);
+
+  DofAffineConstraint<double> cons;
+  for (int i=0; i<_ndofs; i++){
+    for (int j=0; j<_v.size(); j++){
+      if (j!= _positive)
+        cons.linear.push_back(std::pair<Dof,double>(k[j][i],-1.*_weight(j)));
+    };
+    cons.shift=g(i);
+    con[kp[i]] = cons;
+    cons.linear.clear();
+  };
+};
+
+void periodicSupplementConstraint::print() const{
+  Msg::Info("periodicSupplementConstraint Positive = %d",_v[_positive]->getNum());
+  Msg::Info("COnstraint numbered = %d",_tempVertex->getNum());
+  for (int i=0; i<_v.size(); i++){
+    printf("%d \t",_v[i]->getNum());
+  }
+  printf("\n");
+};
+
+void periodicSupplementConstraint::getKinematicalVector(fullVector<double>& m) const{
+  fullVector<double> kinVec;
+	_mbc->getKinematicalVector(kinVec);
+	m.resize(_ndofs);
+	m.setAll(0.);
+	_S.mult(kinVec,m);
+};
+void periodicSupplementConstraint::getKinematicalMatrix(fullMatrix<double>& m) const{
+  m = _S;
+};
+void periodicSupplementConstraint::getLinearConstraints(std::map<Dof,DofAffineConstraint<double> >& con) const{
+
+  fullVector<double> g;
+  getKinematicalVector(g);
+  getLinearConstraints(g, con);
+};
+
+
+
+
 lagrangeSupplementConstraint::lagrangeSupplementConstraint(nonLinearMicroBC* mbc, FunctionSpaceBase* space, FunctionSpaceBase* mspace, 
 												const int c, std::vector<MVertex*>& v): 
 												supplementConstraintBase(mbc,space,mspace,c),

NonLinearSolver/periodicBC/pbcSupplementConstraint.h

@@ -126,6 +126,7 @@ class supplementConstraint : public supplementConstraintBase{
 												const int c,groupOfElements* g, const scalarWeightFunction* func=NULL);
     virtual ~supplementConstraint(){
       if (_weightFunction != NULL)  delete _weightFunction; _weightFunction = NULL;
+      if (_tempVertex) delete _tempVertex;
     }
     virtual void getConstraintKeys(std::vector<Dof>& key) const ; // real dofs on constraint elements
     virtual void getMultiplierKeys(std::vector<Dof>& key) const ;	// multiplier dof on constraint element
@@ -152,6 +153,49 @@ class supplementConstraint : public supplementConstraintBase{
 		virtual void getLinearConstraints(std::map<Dof,DofAffineConstraint<double> >& con) const;
 };
 
+class periodicSupplementConstraint : public supplementConstraintBase{
+	protected:
+    groupOfElements* _gPlus, *_gMinus;
+    fullMatrix<double>  _C, _S, _Cbc; // constraint matrix
+    std::vector<MVertex*> _v; // all Vertex
+    fullVector<double> _weight;
+    int _positive; // positive position
+    SVector3 _Xmean; // int_S{X} dA
+		STensor3 _XXmean; // int_S{0.5*X*X} dA
+    MVertex* _tempVertex;
+    
+	public:
+    periodicSupplementConstraint(nonLinearMicroBC* Mbc, FunctionSpaceBase* space, FunctionSpaceBase* mspace,
+																	const int c, groupOfElements* gplus, groupOfElements* gMinus);
+    virtual ~periodicSupplementConstraint(){
+      if (_tempVertex) delete _tempVertex;
+    }
+    virtual void getConstraintKeys(std::vector<Dof>& key) const ; // real dofs on constraint elements
+    virtual void getMultiplierKeys(std::vector<Dof>& key) const ;	// multiplier dof on constraint element
+    virtual void getConstraintMatrix(fullMatrix<double>& m) const ;		// matrix C
+    virtual void getDependentMatrix(fullMatrix<double>& m) const;
+    virtual void getDependentKeys(std::vector<Dof>& keys) const; // left real dofs
+    virtual void getIndependentKeys(std::vector<Dof>& keys) const; // for spline interpolation
+    virtual void getLinearConstraints(const fullVector<double>& g,
+                   std::map<Dof,DofAffineConstraint<double> >& con) const;
+    virtual constraintElement::ElementType getType() const{return constraintElement::Supplement;};
+    virtual void print() const;
+    
+    virtual std::vector<MVertex*>& getVertexList() {return _v;};
+    virtual fullVector<double>& getWeight() {return _weight;};
+    
+    virtual bool isActive(const MVertex* v) const {
+      for (int i=0; i<_v.size(); i++){
+        if (_v[i]->getNum() == v->getNum()) return true;
+      }
+      return false;
+    };
+
+    virtual void getKinematicalVector(fullVector<double>& m) const;
+		virtual void getKinematicalMatrix(fullMatrix<double>& m) const;
+		virtual void getLinearConstraints(std::map<Dof,DofAffineConstraint<double> >& con) const;
+};
+
 class lagrangeSupplementConstraint : public supplementConstraintBase{
   protected:
     std::vector<MVertex*>& _v; // all Vertex
@@ -169,7 +213,9 @@ class lagrangeSupplementConstraint : public supplementConstraintBase{
   public:
     lagrangeSupplementConstraint(nonLinearMicroBC*, FunctionSpaceBase* space, FunctionSpaceBase* mspace, 
 												const int c, std::vector<MVertex*>& v);
-    virtual ~lagrangeSupplementConstraint(){};
+    virtual ~lagrangeSupplementConstraint(){
+      if (_tempVertex) delete _tempVertex;
+    };
     virtual void getConstraintKeys(std::vector<Dof>& key) const ; // real dofs on constraint elements
     virtual void getMultiplierKeys(std::vector<Dof>& key) const ;	// multiplier dof on constraint element
     virtual void getConstraintMatrix(fullMatrix<double>& m) const ;		// matrix C
@@ -221,7 +267,9 @@ class cubicSplineSupplementConstraint : public supplementConstraintBase{
   public:
     cubicSplineSupplementConstraint(nonLinearMicroBC* mbc, FunctionSpaceBase* space, FunctionSpaceBase* mspace, 
 												const int c, std::vector<MVertex*>& v, int flag = 0);
-    virtual ~cubicSplineSupplementConstraint(){};
+    virtual ~cubicSplineSupplementConstraint(){
+      if (_tempVertex) delete _tempVertex;
+    };
     virtual void getConstraintKeys(std::vector<Dof>& key) const; // real dofs on constraint elements
     virtual void getMultiplierKeys(std::vector<Dof>& key) const;	// multiplier dof on constraint element
     virtual void getConstraintMatrix(fullMatrix<double>& m) const;		// matrix C
@@ -263,7 +311,9 @@ class CoonsPatchLagrangeSupplementConstraint : public supplementConstraintBase{
     CoonsPatchLagrangeSupplementConstraint(nonLinearMicroBC* mbc, FunctionSpaceBase* space, FunctionSpaceBase* mspace,
 																					const int c, lagrangeSupplementConstraint* conX,
                                            lagrangeSupplementConstraint* conY);
-    virtual ~CoonsPatchLagrangeSupplementConstraint(){};
+    virtual ~CoonsPatchLagrangeSupplementConstraint(){
+      if (_tempVertex) delete _tempVertex;
+    };
     virtual void getConstraintKeys(std::vector<Dof>& key) const; // real dofs on constraint elements
     virtual void getMultiplierKeys(std::vector<Dof>& key) const;	// multiplier dof on constraint element
     virtual void getConstraintMatrix(fullMatrix<double>& m) const;		// matrix C