diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
index 4e1e963974ff1024cc565000a31b80b7c90b075d..e82cc641628def330c5216d1afb379c5101c9e6e 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
@@ -262,12 +262,11 @@ nonLinearPeriodicBC::nonLinearPeriodicBC(const int tag, const int dim):
nonLinearMicroBC(tag,dim),_XaddNewVertices(false),_Xdegree(3),
_YaddNewVertices(false),_Ydegree(3),
_ZaddNewVertices(false),_Zdegree(3),
- _wM(nonLinearPeriodicBC::CEM),_forceInterp(false),
- _isShifted(false),_shiftPBCNormal(0.,0.,0.){ }
+ _wM(nonLinearPeriodicBC::CEM),_forceInterp(false){ }
nonLinearPeriodicBC::nonLinearPeriodicBC(const nonLinearPeriodicBC& src):nonLinearMicroBC(src),
_wM(src._wM),_Xdegree(src._Xdegree),_XaddNewVertices(src._XaddNewVertices),_Ydegree(src._Ydegree),
_YaddNewVertices(src._YaddNewVertices),_Zdegree(src._Zdegree),_ZaddNewVertices(src._ZaddNewVertices),
- _forceInterp(src._forceInterp),_isShifted(src._isShifted),_shiftPBCNormal(src._shiftPBCNormal)
+ _forceInterp(src._forceInterp)
{
};
@@ -325,27 +324,50 @@ void nonLinearPeriodicBC::setPeriodicBCOptions(const int method, const int seg,
};
-void nonLinearPeriodicBC::setShiftedBC(const bool fl){
- _isShifted = fl;
-}
-void nonLinearPeriodicBC::setShiftPBCNormal(const double n1, const double n2, const double n3){
+
+
+nonLinearShiftedPeriodicBC::nonLinearShiftedPeriodicBC(const int tag, const int dim): nonLinearPeriodicBC(tag,dim),
+ _shiftPBCNormal(0.,0.,0.){};
+
+nonLinearShiftedPeriodicBC::nonLinearShiftedPeriodicBC(const nonLinearShiftedPeriodicBC& src): nonLinearPeriodicBC(src)
+ ,_shiftPBCNormal(src._shiftPBCNormal){};
+
+
+void nonLinearShiftedPeriodicBC::setShiftPBCNormal(const double n1, const double n2, const double n3){
_shiftPBCNormal[0] = n1;
_shiftPBCNormal[1] = n2;
_shiftPBCNormal[2] = n3;
printf("set Shifted PBC = [%f %f %f]\n",n1,n2,n3);
- _shiftPBCNormal.normalize();
+ if (_shiftPBCNormal.norm()> 0.){
+ _shiftPBCNormal.normalize();
+ }
};
-void nonLinearPeriodicBC::setShiftPBCNormal(const SVector3& n){
+void nonLinearShiftedPeriodicBC::setShiftPBCNormal(const SVector3& n){
this->setShiftPBCNormal(n[0],n[1],n[2]);
printf("set Shifted PBC = [%f %f %f]\n",n[0],n[1],n[2]);
};
-void nonLinearPeriodicBC::updateWithSolver(const nonLinearMechSolver* s){
+void nonLinearShiftedPeriodicBC::updateWithSolver(const nonLinearMechSolver* s){
_shiftPBCNormal = s->getLostSolutionUniquenssNormal();
printf("shitfed with lost of solution uniqueness normal : n = [%f %f %f]",_shiftPBCNormal[0],_shiftPBCNormal[1],_shiftPBCNormal[2]);
};
+nonLinearGeneralPeriodicBC::nonLinearGeneralPeriodicBC(const int tag, const int dim):nonLinearPeriodicBC(tag,dim),
+ _pbcNormal(0.){};
+nonLinearGeneralPeriodicBC::nonLinearGeneralPeriodicBC(const nonLinearGeneralPeriodicBC& src):nonLinearPeriodicBC(src),
+ _pbcNormal(src._pbcNormal){};
+
+void nonLinearGeneralPeriodicBC::setPBCNormal(const double n1, const double n2, const double n3){
+ _pbcNormal[0] = n1;
+ _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):
nonLinearMicroBC(tag,dim){};
@@ -388,58 +410,3 @@ void nonLinearOrthogonalMixedBCDirectionFollowing::setSUBCDirection(const double
void nonLinearOrthogonalMixedBCDirectionFollowing::clearDirs(){ _KUBCDirection.clear(); _SUBCDirection.clear();};
-
-nonLinearMixedBCDirectionFollowing::nonLinearMixedBCDirectionFollowing(const int tag, const int dim):
- nonLinearMicroBC(tag,dim){}
-
-nonLinearMixedBCDirectionFollowing::nonLinearMixedBCDirectionFollowing(const nonLinearMixedBCDirectionFollowing& src):
- nonLinearMicroBC(src),_KUBCDirection(src._KUBCDirection),_SUBCDirection(src._SUBCDirection),
- _PBCNormalDirection(src._PBCNormalDirection), _PBCDirection(src._PBCDirection),
- _rootPhysical(src._rootPhysical),
- _interpolationDegree(src._interpolationDegree),
- _fullConstrained(src._fullConstrained){}
-
-void nonLinearMixedBCDirectionFollowing::setKUBCDirection(const SVector3& dir){
- _KUBCDirection.push_back(dir);
-};
-void nonLinearMixedBCDirectionFollowing::setSUBCDirection(const SVector3& dir){
- _SUBCDirection.push_back(dir);
-};
-
-void nonLinearMixedBCDirectionFollowing::setPBCDirection(const SVector3& dir, const SVector3& norm, const bool full){
- _PBCDirection.push_back(dir);
- _PBCNormalDirection.push_back(norm);
- _fullConstrained.push_back(full);
-};
-
-void nonLinearMixedBCDirectionFollowing::setKUBCDirection(const double nx, const double ny, const double nz){
- SVector3 n(nx,ny,nz);
- n.normalize();
- _KUBCDirection.push_back(n);
-};
-void nonLinearMixedBCDirectionFollowing::setSUBCDirection(const double nx, const double ny, const double nz){
- SVector3 n(nx,ny,nz);
- n.normalize();
- _SUBCDirection.push_back(n);
-};
-
-void nonLinearMixedBCDirectionFollowing::setPBCDirection(const double nx, const double ny, const double nz,
- const double tx, const double ty, const double tz, const bool fullConstrained ){
- SVector3 n(nx,ny,nz);
- n.normalize();
- _PBCNormalDirection.push_back(n);
-
- SVector3 t(tx,ty,tz);
- t.normalize();
- _PBCDirection.push_back(t);
-
- _fullConstrained.push_back(fullConstrained);
-};
-
-void nonLinearMixedBCDirectionFollowing::setRootPhysical(const int phy){
- _rootPhysical.push_back(phy);
-};
-
-void nonLinearMixedBCDirectionFollowing::setInterpolationDegree(const int deg){
- _interpolationDegree.push_back(deg);
-};
diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
index c7cfcc42481405439862b5191eabdf548b669d84..67f607ae7f8e7684ccc907e2ec3be0005d6e31dd 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
@@ -22,7 +22,7 @@ class nonLinearMechSolver;
class nonLinearMicroBC{
public:
- enum whichBCType{LDBC =0, PBC=1, MKBC=2, MIXBC=3, OrthogonalDirectionalMixedBC =4, DirectionalMIXBC=5}; // Linear displacement BC, periodic BC, minimal kinematical BC, mixed BC
+ enum whichBCType{LDBC =0, PBC=1, MKBC=2, MIXBC=3, OrthogonalDirectionalMixedBC =4, GeneralPBC=5, ShiftedPBC = 6}; // Linear displacement BC, periodic BC, minimal kinematical BC, mixed BC
#ifndef SWIG
protected:
@@ -402,7 +402,7 @@ class nonLinearPeriodicBC: public nonLinearMicroBC{
public:
enum whichMethod{CEM=0, LIM=1, CSIM=2, FE_LIN=3, FE_QUA=4, PROJECT=5}; // constraint elimination, lagrange interpolation, cubic spline interpolation --> for constructing linear constraint
- private:
+ protected:
whichMethod _wM;
int _Xdegree; // for lagrange and spline interpolation
bool _XaddNewVertices; // add new nodes for formultation
@@ -417,10 +417,6 @@ class nonLinearPeriodicBC: public nonLinearMicroBC{
bool _forceInterp; // true if forcing setting value
// false if checking with existing vertex
- // for shiftedBC
- bool _isShifted;
- SVector3 _shiftPBCNormal; // by a plane whose normal is _shiftPBCNormal
-
public:
nonLinearPeriodicBC(const int tag, const int dim);
@@ -459,12 +455,6 @@ class nonLinearPeriodicBC: public nonLinearMicroBC{
else if (i==2) return _Zdegree;
else Msg::Fatal("direction must be correctly defined getPBCPolynomialDegree(i)");
}
-
- const SVector3& getShiftPBCNormal() const {return _shiftPBCNormal;};
- void setShiftPBCNormal(const SVector3& n);
- bool isShifted() const{return _isShifted;};
-
- virtual void updateWithSolver(const nonLinearMechSolver* s);
#endif //SWIG
void setAddNewVertices(const int i, const bool add = true) {
@@ -490,10 +480,52 @@ class nonLinearPeriodicBC: public nonLinearMicroBC{
_wM = whichMethod(i);
}
void setPeriodicBCOptions(const int method, const int seg = 5, const bool addvertex = 0);
- void setShiftedBC(const bool fl);
+
+};
+
+class nonLinearShiftedPeriodicBC : public nonLinearPeriodicBC{
+ protected:
+ #ifndef SWIG
+ // for shiftedBC
+ SVector3 _shiftPBCNormal; // by a plane whose normal is _shiftPBCNormal
+ #endif //SWIG
+
+ public:
+ nonLinearShiftedPeriodicBC(const int tag, const int dim);
+ #ifndef SWIG
+ nonLinearShiftedPeriodicBC(const nonLinearShiftedPeriodicBC& src);
+ virtual ~nonLinearShiftedPeriodicBC(){};
+
+ virtual nonLinearMicroBC::whichBCType getType() const {return nonLinearMicroBC::ShiftedPBC;};
+ virtual nonLinearMicroBC* clone() const {return new nonLinearShiftedPeriodicBC(*this);};
+
+ const SVector3& getShiftPBCNormal() const {return _shiftPBCNormal;};
+ void setShiftPBCNormal(const SVector3& n);
+
+ virtual void updateWithSolver(const nonLinearMechSolver* s);
+ #endif
void setShiftPBCNormal(const double n1, const double n2, const double n3);
};
+class nonLinearGeneralPeriodicBC : public nonLinearPeriodicBC{
+ protected:
+ #ifndef SWIG
+ SVector3 _pbcNormal; // unit normal of PBC plane
+ #endif //SWIG
+ public:
+ nonLinearGeneralPeriodicBC(const int tag, const int dim);
+ #ifndef SWIG
+ nonLinearGeneralPeriodicBC(const nonLinearGeneralPeriodicBC& src);
+ virtual ~nonLinearGeneralPeriodicBC(){};
+ virtual nonLinearMicroBC::whichBCType getType() const {return nonLinearMicroBC::GeneralPBC;};
+ virtual nonLinearMicroBC* clone() const {return new nonLinearGeneralPeriodicBC(*this);};
+
+ const SVector3& getPBCNormal() const {return _pbcNormal;};
+ void setPBCNormal(const SVector3& n) {_pbcNormal = n;}
+ #endif //SWIG
+ void setPBCNormal(const double n1, const double n2, const double n3);
+};
+
class nonLinearDisplacementBC : public nonLinearMicroBC{
@@ -613,53 +645,5 @@ class nonLinearOrthogonalMixedBCDirectionFollowing : public nonLinearMicroBC{
void setSUBCDirection(const double nx, const double ny, const double nz);
};
-class nonLinearMixedBCDirectionFollowing : public nonLinearMicroBC{
- #ifndef SWIG
- protected:
- std::vector<SVector3> _KUBCDirection;
- std::vector<SVector3> _SUBCDirection;
- std::vector<SVector3> _PBCDirection;
- std::vector<SVector3> _PBCNormalDirection; // normal to PBCdir
-
- std::vector<int> _rootPhysical; // v1 v2 v4 for 2D and v1 v2 v4 v5 for 3D
- std::vector<int> _interpolationDegree; // following v1v2, v1v4, v1v5
- std::vector<bool> _fullConstrained; // true if full, false otherwise
-
- #endif // SWIG
- public:
- nonLinearMixedBCDirectionFollowing(const int tag, const int dim);
- #ifndef SWIG
- nonLinearMixedBCDirectionFollowing(const nonLinearMixedBCDirectionFollowing& src);
- virtual ~nonLinearMixedBCDirectionFollowing(){};
- virtual nonLinearMicroBC::whichBCType getType() const{
- return nonLinearMicroBC::DirectionalMIXBC;
- };
-
- virtual nonLinearMicroBC* clone() const {return new nonLinearMixedBCDirectionFollowing(*this);}
-
- const std::vector<SVector3>& getKUBCDirection() const {return _KUBCDirection;};
- const std::vector<SVector3>& getSUBCDirection() const {return _SUBCDirection;};
-
- // for PBC
- const std::vector<SVector3>& getPBCNormalDirection() const {return _PBCNormalDirection;};
- const std::vector<SVector3>& getPBCDirection() const {return _PBCDirection;};
- const std::vector<int>& getRootPhysical() const {return _rootPhysical;};
- const std::vector<int>& getInterpolationDegree() const {return _interpolationDegree;};
- const std::vector<bool>& getFullConstrainedFlag() const {return _fullConstrained;};
-
- void setKUBCDirection(const SVector3& dir);
- void setSUBCDirection(const SVector3& dir);
- void setPBCDirection(const SVector3& dir, const SVector3& norm, const bool full);
-
- #endif //SWIG
-
- void setKUBCDirection(const double nx, const double ny, const double nz);
- void setSUBCDirection(const double nx, const double ny, const double nz);
- void setPBCDirection(const double nx, const double ny, const double nz,
- const double tx, const double ty, const double tz,
- const bool fullConstrained);
- void setRootPhysical(const int phy);
- void setInterpolationDegree(const int deg);
-};
#endif // NONLINEARMICROBC_H_
diff --git a/NonLinearSolver/periodicBC/directionalConstraintElement.cpp b/NonLinearSolver/periodicBC/directionalConstraintElement.cpp
index 8ecaf5d5b65ce018b99d3c8060717a0f2ab3f101..16c7162bc0eaf82c4f24a0c6411e8982da8b05b6 100644
--- a/NonLinearSolver/periodicBC/directionalConstraintElement.cpp
+++ b/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),
diff --git a/NonLinearSolver/periodicBC/pbcConstraintElement.cpp b/NonLinearSolver/periodicBC/pbcConstraintElement.cpp
index f84b5cc020cd10fe983aad60245142654730f8ec..a11d8751586c817b6a6ef4afb97ff01e64bb1bc7 100644
--- a/NonLinearSolver/periodicBC/pbcConstraintElement.cpp
+++ b/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]);
};
diff --git a/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp b/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
index 816d95a90707485c8eeee8e00d3d6411f47b1bb1..60778c8dcdf98623adb1f2cba862d30938ae475c 100644
--- a/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
+++ b/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
@@ -1676,16 +1676,9 @@ void pbcConstraintElementGroup::createCubicSplineConstraintElementGroup(){
};
void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
- nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
+ nonLinearShiftedPeriodicBC* pbc = static_cast<nonLinearShiftedPeriodicBC*>(_solver->getMicroBC());
const SVector3& normal = pbc->getShiftPBCNormal();
- normal.print("normal");
- if (normal.norm() < 1e-10){
- printf("no shifted if normal is zero pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup \n");
- this->createLagrangeConstraintElementGroup();
- pbc->setShiftedBC(false);
- return;
- }
int dim = _solver->getMicroBC()->getDim();
const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
@@ -1735,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();
@@ -1803,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(){
@@ -2004,123 +2070,7 @@ void pbcConstraintElementGroup::createOrthogonalMixBCConstraintElementGroupDirec
}
};
-void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollowing(){
- nonLinearMixedBCDirectionFollowing* mixBC = dynamic_cast<nonLinearMixedBCDirectionFollowing*>(_solver->getMicroBC());
- const std::vector<SVector3>& kubcDir = mixBC->getKUBCDirection();
- const std::vector<SVector3>& subcDir = mixBC->getSUBCDirection();
- const std::vector<SVector3>& pbcDir = mixBC->getPBCDirection();
- const std::vector<SVector3>& pbcNormalDir = mixBC->getPBCNormalDirection();
- const std::vector<bool>& pbcFullConstrainedFlag = mixBC->getFullConstrainedFlag();
-
-
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- std::set<MVertex*> allVertex;
-
- for (int i=0; i< bgroup.size(); i++){
- for (groupOfElements::vertexContainer::iterator it = bgroup[i]->vbegin(); it!= bgroup[i]->vend(); it++){
- allVertex.insert(*it);
- }
- }
-
- for (int i=0; i< kubcDir.size(); i++){
- const SVector3& dir = kubcDir[i];
- for (std::set<MVertex*>::iterator it = allVertex.begin(); it!= allVertex.end(); it++){
- MVertex* v = *it;
-
- if ((pbcNormalDir.size() == 0) and (it == allVertex.begin())){
- MVertex* v = *(allVertex.begin());
- constraintElement* ppc = new periodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,-1,v);
- _allConstraint.insert(ppc);
- }
- else{
- constraintElement* ppc = new directionalPeriodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,v,NULL,dir);
- _allConstraint.insert(ppc);
- }
- }
- }
- if (pbcNormalDir.size() > 0){
- const std::vector<int>& rootPhysical = mixBC->getRootPhysical();
- const std::vector<int>& interpDeg = mixBC->getInterpolationDegree();
-
- if (mixBC->getDim() == 2){
- MVertex* v1(NULL), *v2(NULL), *v4(NULL);
- for (std::set<MVertex*>::iterator it = allVertex.begin(); it!= allVertex.end(); it++){
- MVertex* v = *it;
- if ((v1 == NULL) and (v->getNum() == rootPhysical[0])){
- v1 = v;
- }
- if ((v2 == NULL) and (v->getNum() == rootPhysical[1])){
- v2 = v;
- }
- if ((v4 == NULL) and (v->getNum() == rootPhysical[2])){
- v4 = v;
- }
- if ((v1 != NULL) and (v2!=NULL) and (v4!=NULL)){
- break;
- }
- };
- if ((v1==NULL) or (v2 == NULL) or (v4 ==NULL)){
- Msg::Fatal("root physical vertces are not correctly defined");
- }
-
- constraintElement* ppc = new periodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,-1,v2);
- _allConstraint.insert(ppc);
-
- // create local basi
- InterpolationOperations::createLinePolynomialBase(interpDeg[0],false,_v1,_v2,l12,_xBase,_virtualVertices);
- InterpolationOperations::createLinePolynomialBase(interpDeg[0],false,_v1,_v4,l41,_yBase,_virtualVertices);
-
- for (int i=0; i< _xBase.size(); i++){
- Msg::Info("xbase : %d",_xBase[i]->getNum());
- }
-
- for (int i=0; i< _yBase.size(); i++){
- Msg::Info("ybase : %d",_yBase[i]->getNum());
- }
-
- std::set<MVertex*> allPositiveVertex;
- for (groupOfElements::vertexContainer::iterator it = bgroup[2]->vbegin(); it!= bgroup[2]->vend(); it++){
- allPositiveVertex.insert(*it);
- }
- for (groupOfElements::vertexContainer::iterator it = bgroup[3]->vbegin(); it!= bgroup[3]->vend(); it++){
- allPositiveVertex.insert(*it);
- }
-
-
- for (int i=0; i< pbcNormalDir.size(); i++){
- const SVector3& normalPBC = pbcNormalDir[i];
- const SVector3& dirPBC = pbcDir[i];
- const bool& fullCon = pbcFullConstrainedFlag[i];
-
- if (fullCon){
- this->lagrangeFormCondition(bgroup[0]->vbegin(),bgroup[0]->vend(),_xBase);
- this->lagrangeFormCondition(bgroup[1]->vbegin(),bgroup[1]->vend(),_yBase);
- this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,normalPBC,true);
- }
- else{
- this->directionalLagrangeFormCondition(bgroup[0]->vbegin(),bgroup[0]->vend(),_xBase,dirPBC);
- this->directionalLagrangeFormCondition(bgroup[1]->vbegin(),bgroup[1]->vend(),_yBase,dirPBC);
- this->directionalLagrangePeriodicCondition(allPositiveVertex.begin(), allPositiveVertex.end(),_xBase,_yBase,normalPBC,false);
- }
- }
- }
- else{
- Msg::Fatal("this bc is not implemented for 3D problems");
- }
- }
-
- int halfOfgsize = bgroup.size()/2;
- for (int i=0; i< subcDir.size(); i++){
- const SVector3& dir = subcDir[i];
- for (int j=0; j< halfOfgsize; j++){
- groupOfElements* gPlus = bgroup[j];
- groupOfElements* gMinus = bgroup[j+halfOfgsize];
- constraintElement* ppc = new directionalPBCSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,gPlus,gMinus,dir);
- _allConstraint.insert(ppc);
- };
- }
-};
void pbcConstraintElementGroup::createMixBCConstraintElementGroup_DG(){
nonLinearMixedBC* mixBC = dynamic_cast<nonLinearMixedBC*>(_solver->getMicroBC());
@@ -2219,397 +2169,508 @@ void pbcConstraintElementGroup::createCornerConstraintElementGroupShiftedBPC(){
if (_solver->getMicroBC()->getOrder() == 2){
Msg::Fatal("shifted BC is implemented for first order BC only");
}
+
+ nonLinearShiftedPeriodicBC* spbc = static_cast<nonLinearShiftedPeriodicBC*>(_solver->getMicroBC());
- if (_cVertices.size() > 0){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,-1,_v1);
- _allConstraint.insert(cel);
- }
+ if (spbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
+ if (_cVertices.size() > 0){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,-1,_v1);
+ _allConstraint.insert(cel);
+ }
+ }
+ else{
+ Msg::Fatal("this method is not implemented to impose shiftedPBC");
+ }
+
+
}
void pbcConstraintElementGroup::createCornerConstraintElementGroup(){
-
- if (_solver->getMicroBC()->getType() == nonLinearMicroBC::PBC){
- int ndofsPerNode = _solver->getMicroBC()->getNumberOfConstrainedDofsPerVertex();
- std::vector<int> mechanics(_solver->getMicroBC()->getMechanicsConstrainedComps().begin(), _solver->getMicroBC()->getMechanicsConstrainedComps().end()); // mechanic Dof only
- std::vector<int> compConstitutive(_solver->getMicroBC()->getConstitutiveExtraConstrainedComps().begin(), _solver->getMicroBC()->getConstitutiveExtraConstrainedComps().end()); // constitutuive extraDofs
- std::vector<int> compNotConstitutive(_solver->getMicroBC()->getNonConstitutiveExtraConstrainedComps().begin(), _solver->getMicroBC()->getNonConstitutiveExtraConstrainedComps().end()); // non-constitutive extraDofs
+ nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
+ int ndofsPerNode = _solver->getMicroBC()->getNumberOfConstrainedDofsPerVertex();
+ std::vector<int> mechanics(_solver->getMicroBC()->getMechanicsConstrainedComps().begin(), _solver->getMicroBC()->getMechanicsConstrainedComps().end()); // mechanic Dof only
+ std::vector<int> compConstitutive(_solver->getMicroBC()->getConstitutiveExtraConstrainedComps().begin(), _solver->getMicroBC()->getConstitutiveExtraConstrainedComps().end()); // constitutuive extraDofs
+ std::vector<int> compNotConstitutive(_solver->getMicroBC()->getNonConstitutiveExtraConstrainedComps().begin(), _solver->getMicroBC()->getNonConstitutiveExtraConstrainedComps().end()); // non-constitutive extraDofs
- nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::CEM or pbc->getPBCMethod() == nonLinearPeriodicBC::PROJECT){
- // for mechanical Dofs
- std::set<MVertex*> activeCornerVertices;
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
- if (!isVirtual(vp1) or (isVirtual(vp1) and isActive(vp1))){
- activeCornerVertices.insert(vp1);
- }
- }
- if (activeCornerVertices.size() == 0){
- Msg::Info("all corner vertices are not active");
+ const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
+ if (pbc->getPBCMethod() == nonLinearPeriodicBC::CEM or pbc->getPBCMethod() == nonLinearPeriodicBC::PROJECT){
+ // for mechanical Dofs
+ std::set<MVertex*> activeCornerVertices;
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ if (!isVirtual(vp1) or (isVirtual(vp1) and isActive(vp1))){
+ activeCornerVertices.insert(vp1);
}
-
- if (_solver->getMicroBC()->getOrder() == 1){
- // fix all corners mechanical Dofs
- if (activeCornerVertices.size() > 0){
- for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
- MVertex* vp1 = *it;
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
- _allConstraint.insert(cel);
- }
- }
- }
- else{
- // fix an arbitrary point
- MVertex* v = *(bgroup[0]->vbegin());
+ }
+ if (activeCornerVertices.size() == 0){
+ Msg::Info("all corner vertices are not active");
+ }
+
+ if (_solver->getMicroBC()->getOrder() == 1){
+ // fix all corners mechanical Dofs
+ if (activeCornerVertices.size() > 0){
+ for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ MVertex* vp1 = *it;
for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v);
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
_allConstraint.insert(cel);
}
- }
+ }
}
- else if (_solver->getMicroBC()->getOrder() == 2){
- // periodic between
- if (activeCornerVertices.size() >0){
- MVertex* vn = *(activeCornerVertices.begin());
- for(std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
- if (it != activeCornerVertices.begin())
- {
- MVertex* vp1 = *it;
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i], vp1,vn);
- _allConstraint.insert(cel);
- }
- }
- }
- }
-
- for (int i=0; i<bgroup.size()/2; i++){
- groupOfElements* g = bgroup[i];
- for (int k=0; k< mechanics.size(); k++){
- constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],g);
- _allConstraint.insert(ppc);
- }
+ else{
+ // fix an arbitrary point
+ MVertex* v = *(bgroup[0]->vbegin());
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v);
+ _allConstraint.insert(cel);
}
}
-
- // for constitutive extraDofs
- if (activeCornerVertices.size() > 0){
- MVertex* oneV = *(activeCornerVertices.begin());
- for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
- MVertex* vp1 = *it;
- if (vp1 != oneV){
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,oneV);
+ }
+ else if (_solver->getMicroBC()->getOrder() == 2){
+ // periodic between
+ if (activeCornerVertices.size() >0){
+ MVertex* vn = *(activeCornerVertices.begin());
+ for(std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ if (it != activeCornerVertices.begin())
+ {
+ MVertex* vp1 = *it;
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i], vp1,vn);
_allConstraint.insert(cel);
}
}
}
}
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,compConstitutive[i],_solver->getDomainVector());
- _allConstraint.insert(celExtra);
+
+ for (int i=0; i<bgroup.size()/2; i++){
+ groupOfElements* g = bgroup[i];
+ for (int k=0; k< mechanics.size(); k++){
+ constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],g);
+ _allConstraint.insert(ppc);
+ }
}
+ }
- // for non-constitutive extraDofs
- if (activeCornerVertices.size() >0){
- for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
- MVertex* vp1 = *it;
- for (int i=0; i< compNotConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
+ // for constitutive extraDofs
+ if (activeCornerVertices.size() > 0){
+ MVertex* oneV = *(activeCornerVertices.begin());
+ for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ MVertex* vp1 = *it;
+ if (vp1 != oneV){
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,oneV);
_allConstraint.insert(cel);
}
- }
+ }
}
- else{
- MVertex* v = *(bgroup[0]->vbegin());
+ }
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,compConstitutive[i],_solver->getDomainVector());
+ _allConstraint.insert(celExtra);
+ }
+
+ // for non-constitutive extraDofs
+ if (activeCornerVertices.size() >0){
+ for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ MVertex* vp1 = *it;
for (int i=0; i< compNotConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],v);
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
_allConstraint.insert(cel);
}
-
+ }
+ }
+ else{
+ MVertex* v = *(bgroup[0]->vbegin());
+ for (int i=0; i< compNotConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],v);
+ _allConstraint.insert(cel);
}
+
}
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM){
- // for mechanics
- if (_solver->getMicroBC()->getOrder() == 1){
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM){
+ // for mechanics
+ if (_solver->getMicroBC()->getOrder() == 1){
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
+ _allConstraint.insert(cel);
+ }
+ }
+ }
+ else if (_solver->getMicroBC()->getOrder() == 2){
+ for(std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* v = *it;
+ if ( v != _v1){
for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v,_v1);
_allConstraint.insert(cel);
}
}
}
- else if (_solver->getMicroBC()->getOrder() == 2){
- for(std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* v = *it;
- if ( v != _v1){
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v,_v1);
- _allConstraint.insert(cel);
- }
- }
- }
- if (_solver->getMicroBC()->getDim() ==2){
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* ppc = new cubicSplineSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_xBase,0);
- _allConstraint.insert(ppc);
+ if (_solver->getMicroBC()->getDim() ==2){
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* ppc = new cubicSplineSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_xBase,0);
+ _allConstraint.insert(ppc);
- ppc = new cubicSplineSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_yBase,1);
+ ppc = new cubicSplineSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_yBase,1);
+ _allConstraint.insert(ppc);
+ }
+ }
+ else if (_solver->getMicroBC()->getDim() == 3){
+ for (int i=0; i<bgroup.size()/2; i++){
+ groupOfElements* g = bgroup[i];
+ for (int k=0; k< mechanics.size(); k++){
+ constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],g);
_allConstraint.insert(ppc);
}
}
- else if (_solver->getMicroBC()->getDim() == 3){
- for (int i=0; i<bgroup.size()/2; i++){
- groupOfElements* g = bgroup[i];
- for (int k=0; k< mechanics.size(); k++){
- constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],g);
- _allConstraint.insert(ppc);
- }
- }
- };
};
+ };
- //
- // for constitutive extraDofs
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
- if (vp1 != _v1){
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,_v1);
- _allConstraint.insert(cel);
- }
+ //
+ // for constitutive extraDofs
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ if (vp1 != _v1){
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,_v1);
+ _allConstraint.insert(cel);
}
}
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
- compConstitutive[i], _solver->getDomainVector());
- _allConstraint.insert(celExtra);
+ }
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
+ compConstitutive[i], _solver->getDomainVector());
+ _allConstraint.insert(celExtra);
+ }
+
+ // for non-constitutive extraDofs
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ for (int i=0; i< compNotConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
+ _allConstraint.insert(cel);
}
- // for non-constitutive extraDofs
+ }
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
+ // for mechanics
+ if (_solver->getMicroBC()->getOrder() == 1){
for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
MVertex* vp1 = *it;
- for (int i=0; i< compNotConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
_allConstraint.insert(cel);
}
-
}
}
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
- // for mechanics
- if (_solver->getMicroBC()->getOrder() == 1){
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
- _allConstraint.insert(cel);
- }
- }
- }
- else if (_solver->getMicroBC()->getOrder() == 2){
- if (pbc->getMaxDegree()>1){
- for(std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* v = *it;
- if ( v != _v1){
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v,_v1);
- _allConstraint.insert(cel);
- }
- }
- }
-
- //
- for (int k=0; k< mechanics.size(); k++){
- if (_solver->getMicroBC()->getDim() ==2){
- constraintElement* ppc = new lagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],_xBase);
- _allConstraint.insert(ppc);
- ppc = new lagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],_yBase);
- _allConstraint.insert(ppc);
- }
- else if (_solver->getMicroBC()->getDim() ==3){
- std::vector<lagrangeSupplementConstraint*> lconstraint;
- lagrangeSupplementConstraint* con = new lagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],_xBase);
- lconstraint.push_back(con);
- con = new lagrangeSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],_yBase);
- lconstraint.push_back(con);
- con = new lagrangeSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],_zBase);
- lconstraint.push_back(con);
- //
- constraintElement* ppc = new CoonsPatchLagrangeSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],lconstraint[0],lconstraint[1]);
- _allConstraint.insert(ppc);
- ppc = new CoonsPatchLagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],lconstraint[1],lconstraint[2]);
- _allConstraint.insert(ppc);
- ppc = new CoonsPatchLagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],lconstraint[0],lconstraint[2]);
- _allConstraint.insert(ppc);
- }
-
- }
- }
- else if (pbc->getMaxDegree() == 1) {
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
+ else if (_solver->getMicroBC()->getOrder() == 2){
+ if (pbc->getMaxDegree()>1){
+ for(std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* v = *it;
+ if ( v != _v1){
for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],vp1);
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v,_v1);
_allConstraint.insert(cel);
}
}
}
- }
+ //
+ for (int k=0; k< mechanics.size(); k++){
+ if (_solver->getMicroBC()->getDim() ==2){
+ constraintElement* ppc = new lagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],_xBase);
+ _allConstraint.insert(ppc);
+ ppc = new lagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],_yBase);
+ _allConstraint.insert(ppc);
+ }
+ else if (_solver->getMicroBC()->getDim() ==3){
+ std::vector<lagrangeSupplementConstraint*> lconstraint;
+ lagrangeSupplementConstraint* con = new lagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],_xBase);
+ lconstraint.push_back(con);
+ con = new lagrangeSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],_yBase);
+ lconstraint.push_back(con);
+ con = new lagrangeSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],_zBase);
+ lconstraint.push_back(con);
+ //
+ constraintElement* ppc = new CoonsPatchLagrangeSupplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[k],lconstraint[0],lconstraint[1]);
+ _allConstraint.insert(ppc);
+ ppc = new CoonsPatchLagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],lconstraint[1],lconstraint[2]);
+ _allConstraint.insert(ppc);
+ ppc = new CoonsPatchLagrangeSupplementConstraint(_solver->getMicroBC(), _LagSpace,_MultSpace,mechanics[k],lconstraint[0],lconstraint[2]);
+ _allConstraint.insert(ppc);
+ }
- // for constitutive extraDofs
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
- if (vp1 != _v1){
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,_v1);
+ }
+ }
+ else if (pbc->getMaxDegree() == 1) {
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],vp1);
_allConstraint.insert(cel);
}
}
}
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
- compConstitutive[i], _solver->getDomainVector());
- _allConstraint.insert(celExtra);
- }
- // for non-constitutive extraDofs
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
- for (int i=0; i< compNotConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
+ }
+
+ // for constitutive extraDofs
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ if (vp1 != _v1){
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,_v1);
_allConstraint.insert(cel);
}
}
+ }
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
+ compConstitutive[i], _solver->getDomainVector());
+ _allConstraint.insert(celExtra);
+ }
+ // for non-constitutive extraDofs
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ for (int i=0; i< compNotConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
+ _allConstraint.insert(cel);
+ }
}
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN or pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
- std::set<MVertex*> activeCornerVertices;
- // all vertual vertices and not appears in any constraint element do not neen to take into account
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* vp1 = *it;
- if (!isVirtual(vp1) or (isVirtual(vp1) and isActive(vp1))){
- activeCornerVertices.insert(vp1);
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN or pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
+
+ std::set<MVertex*> activeCornerVertices;
+ // all vertual vertices and not appears in any constraint element do not neen to take into account
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ MVertex* vp1 = *it;
+ if (!isVirtual(vp1) or (isVirtual(vp1) and isActive(vp1))){
+ activeCornerVertices.insert(vp1);
+ }
+ }
+
+ MVertex* vertexToFix = NULL;
+
+ if (activeCornerVertices.size() == 0){
+ // need to fix some vertices to avoid rigid body mode
+ Msg::Info("all corner vertices are virtual and not active");
+ for (groupOfElements::vertexContainer::iterator itg = _gX->vbegin(); itg!= _gX->vend(); itg++){
+ MVertex* vg = *itg;
+ if (this->isActive(vg)){
+ vertexToFix = vg;
+ break;
}
}
-
- MVertex* vertexToFix = NULL;
-
- if (activeCornerVertices.size() == 0){
- // need to fix some vertices to avoid rigid body mode
- Msg::Info("all corner vertices are virtual and not active");
- for (groupOfElements::vertexContainer::iterator itg = _gX->vbegin(); itg!= _gX->vend(); itg++){
- MVertex* vg = *itg;
- if (this->isActive(vg)){
- vertexToFix = vg;
- break;
+ }
+
+ // for mechanics
+ if (_solver->getMicroBC()->getOrder() == 1){
+ if (activeCornerVertices.size() > 0){
+ for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ MVertex* vp1 = *it;
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
+ _allConstraint.insert(cel);
}
+ }
+ }
+ else{
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vertexToFix);
+ _allConstraint.insert(cel);
}
}
- // for mechanics
- if (_solver->getMicroBC()->getOrder() == 1){
- if (activeCornerVertices.size() > 0){
- for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ }
+ else if (_solver->getMicroBC()->getOrder() == 2){
+ if (activeCornerVertices.size() >0){
+ MVertex* vn = *(activeCornerVertices.begin());
+ for(std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ if (it != activeCornerVertices.begin())
+ {
MVertex* vp1 = *it;
for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vp1);
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],vp1,vn);
_allConstraint.insert(cel);
+
}
- }
- }
- else{
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],vertexToFix);
- _allConstraint.insert(cel);
}
}
-
}
- else if (_solver->getMicroBC()->getOrder() == 2){
- if (activeCornerVertices.size() >0){
- MVertex* vn = *(activeCornerVertices.begin());
- for(std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
- if (it != activeCornerVertices.begin())
- {
- MVertex* vp1 = *it;
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],vp1,vn);
- _allConstraint.insert(cel);
-
- }
- }
- }
+
+ for (int i=0; i< mechanics.size(); i++){
+ if (_solver->getMicroBC()->getDim()== 2){
+ constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gX);
+ _allConstraint.insert(ppc);
+ ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gY);
+ _allConstraint.insert(ppc);
+ }
+ else if (_solver->getMicroBC()->getDim() == 3){
+ constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gXY);
+ _allConstraint.insert(ppc);
+ ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gYZ);
+ _allConstraint.insert(ppc);
+ ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gZX);
+ _allConstraint.insert(ppc);
}
+ }
+ }
- for (int i=0; i< mechanics.size(); i++){
- if (_solver->getMicroBC()->getDim()== 2){
- constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gX);
- _allConstraint.insert(ppc);
- ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gY);
- _allConstraint.insert(ppc);
- }
- else if (_solver->getMicroBC()->getDim() == 3){
- constraintElement* ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gXY);
- _allConstraint.insert(ppc);
- ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gYZ);
- _allConstraint.insert(ppc);
- ppc = new supplementConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],_gZX);
- _allConstraint.insert(ppc);
+ // for constitutive extraDofs
+ if (activeCornerVertices.size() > 0){
+ MVertex* vn1 = *(activeCornerVertices.begin());
+ for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ MVertex* vp1 = *it;
+ if (vp1 != vn1){
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,vn1);
+ _allConstraint.insert(cel);
}
}
}
+ }
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
+ compConstitutive[i], _solver->getDomainVector());
+ _allConstraint.insert(celExtra);
+ }
- // for constitutive extraDofs
- if (activeCornerVertices.size() > 0){
- MVertex* vn1 = *(activeCornerVertices.begin());
+ // for non-constitutive extraDofs
+ if (activeCornerVertices.size() > 0){
+ for (int i=0; i< compNotConstitutive.size(); i++){
for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
MVertex* vp1 = *it;
- if (vp1 != vn1){
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,compConstitutive[i],vp1,vn1);
- _allConstraint.insert(cel);
- }
- }
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
+ _allConstraint.insert(cel);
}
- }
- for (int i=0; i< compConstitutive.size(); i++){
- constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
- compConstitutive[i], _solver->getDomainVector());
- _allConstraint.insert(celExtra);
- }
- // for non-constitutive extraDofs
- if (activeCornerVertices.size() > 0){
- for (int i=0; i< compNotConstitutive.size(); i++){
- for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
- MVertex* vp1 = *it;
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vp1);
- _allConstraint.insert(cel);
- }
+ }
+ }
+ else{
+ for (int i=0; i< compNotConstitutive.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vertexToFix);
+ _allConstraint.insert(cel);
- }
- }
- else{
- for (int i=0; i< compNotConstitutive.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,compNotConstitutive[i],vertexToFix);
- _allConstraint.insert(cel);
+ }
+ }
- }
- }
+ }
+ else
+ Msg::Error("pbcConstraintElementGroup::createCornerConstraintElementGroup: case undefined");
+
+};
+void pbcConstraintElementGroup::createPBCConstraintGroup(){
+ nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
+
+ if (pbc->getPBCMethod() == nonLinearPeriodicBC::CEM){
+ Msg::Info("Imposing PBC by periodic mesh formulation");
+ this->createPeriodicMeshConstraintElementGroup();
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM){
+ Msg::Info("Imposing PBC by cubic spline formulation");
+ this->createCubicSplineConstraintElementGroup();
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
+ Msg::Info("Imposing PBC by lagrange formulation");
+ if (pbc->getMaxDegree() >1){
+ this->createLagrangeConstraintElementGroup();
+ }
+ else{
+ createLinearDisplacementConstraintElementGroup();
}
- else
- Msg::Error("pbcConstraintElementGroup::createCornerConstraintElementGroup: case undefined");
}
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN or pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
+ Msg::Info("Imposing PBC by lagrange C0");
+ this->createFEConstraintElementGroup();
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::PROJECT){
+ Msg::Info("Imposing PBC by projection method");
+ this->createProjectConstraintElementGroup();
+ }
+ else{
+ Msg::Fatal("this method is not implemented to impose PBC");
+ }
+ 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();
+ if (normal.norm() < 1e-10){
+ printf("because shifted vector is null, usual PBC is used\n");
+ this->createPBCConstraintGroup();
+ }
+ else{
+ //
+ printf("PBC is shifted using Lagrange interpolation method\n");
+ spbc->setPBCMethod(nonLinearPeriodicBC::LIM);
+
+ if (spbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
+ Msg::Info("Imposing PBC by lagrange formulation");
+ this->createShiftedLagrangeConstraintElementGroup();
+ }
+ else{
+ Msg::Fatal("this method is not implemented to impose shiftedPBC");
+ };
+
+ this->createCornerConstraintElementGroupShiftedBPC();
+ }
};
@@ -2635,54 +2696,14 @@ void pbcConstraintElementGroup::createConstraintGroup(){
}
else{
if (_solver->getMicroBC()->getType() == nonLinearMicroBC::PBC){
- nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- if (pbc->isShifted() and pbc->getPBCMethod() != nonLinearPeriodicBC::LIM){
- Msg::Warning("shifted has not implemented for this option");
- }
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- if (bgroup.size()>0){
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::CEM){
- Msg::Info("Imposing PBC by periodic mesh formulation");
- this->createPeriodicMeshConstraintElementGroup();
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM){
- Msg::Info("Imposing PBC by cubic spline formulation");
- this->createCubicSplineConstraintElementGroup();
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::LIM){
- Msg::Info("Imposing PBC by lagrange formulation");
- if (pbc->getMaxDegree() >1){
- if (pbc->isShifted()){
- this->createShiftedLagrangeConstraintElementGroup();
- }
- else{
- this->createLagrangeConstraintElementGroup();
- }
- }
- else{
- createLinearDisplacementConstraintElementGroup();
- }
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN or pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
- Msg::Info("Imposing PBC by lagrange C0");
- this->createFEConstraintElementGroup();
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::PROJECT){
- Msg::Info("Imposing PBC by projection method");
- this->createProjectConstraintElementGroup();
- }
- else{
- Msg::Fatal("this method is not implemented to impose PBC");
- }
-
- if (pbc->isShifted()){
- this->createCornerConstraintElementGroupShiftedBPC();
- }
- else{
- this->createCornerConstraintElementGroup();
- }
- }
+ this->createPBCConstraintGroup();
}
+ 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");
this->createLinearDisplacementConstraintElementGroup();
@@ -2695,9 +2716,6 @@ void pbcConstraintElementGroup::createConstraintGroup(){
Msg::Info("imposing mixed BC");
this->createMixBCConstraintElementGroup();
}
- else if (_solver->getMicroBC()->getType() == nonLinearMicroBC::DirectionalMIXBC){
- this->createMixBCConstraintElementGroupDirectionFollowing();
- }
else if (_solver->getMicroBC()->getType() == nonLinearMicroBC::OrthogonalDirectionalMixedBC){
this->createOrthogonalMixBCConstraintElementGroupDirectionFollowing();
}
diff --git a/NonLinearSolver/periodicBC/pbcCreateConstraints.h b/NonLinearSolver/periodicBC/pbcCreateConstraints.h
index 828895456d82e48bc8b2302adcc0ce13fa9ba1aa..23ea069cd8e5b05083ffd52aab12c9d69ffb1922 100644
--- a/NonLinearSolver/periodicBC/pbcCreateConstraints.h
+++ b/NonLinearSolver/periodicBC/pbcCreateConstraints.h
@@ -240,6 +240,7 @@ class pbcConstraintElementGroup{
void createCubicSplineConstraintElementGroup();
void createLagrangeConstraintElementGroup();
void createShiftedLagrangeConstraintElementGroup();
+ void createGeneralPBCLagrangeConstraintElementGroup();
void createFEConstraintElementGroup();
void createLinearDisplacementConstraintElementGroup();
void createProjectConstraintElementGroup();
@@ -255,8 +256,10 @@ class pbcConstraintElementGroup{
void createMixBCConstraintElementGroup_DG();
void createOrthogonalMixBCConstraintElementGroupDirectionFollowing();
- void createMixBCConstraintElementGroupDirectionFollowing();
-
+
+ void createPBCConstraintGroup();
+ void createShiftedPBCConstraintGroup();
+ void createGeneralPBCConstraintGroup();
SVector3 getUniformDisp(MVertex* v);
SVector3 getPerturbation(dofManager<double>* pmanager, MVertex* v);
diff --git a/NonLinearSolver/periodicBC/pbcSupplementConstraint.cpp b/NonLinearSolver/periodicBC/pbcSupplementConstraint.cpp
index 1cec495b39f0d0916bca3f286adf7949db440710..73cb84cadc71d8e0dba71e2d50285b0478ef2012 100644
--- a/NonLinearSolver/periodicBC/pbcSupplementConstraint.cpp
+++ b/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),
diff --git a/NonLinearSolver/periodicBC/pbcSupplementConstraint.h b/NonLinearSolver/periodicBC/pbcSupplementConstraint.h
index 19838a67597473f87fc619256368e1e57de46e1b..38c7cbb0deca9832bd8e18e405d16ed4a1c7720e 100644
--- a/NonLinearSolver/periodicBC/pbcSupplementConstraint.h
+++ b/NonLinearSolver/periodicBC/pbcSupplementConstraint.h
@@ -110,7 +110,6 @@ class supplementConstraint : public supplementConstraintBase{
protected:
groupOfElements* _g;
- QuadratureBase* _integrationRule; // integration rule for approximate integrals
fullMatrix<double> _C, _S, _Cbc; // constraint matrix
SVector3 _Xmean; // int_S{X} dA
STensor3 _XXmean; // int_S{0.5*X*X} dA
@@ -126,6 +125,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 +152,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 +212,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 +266,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 +310,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
diff --git a/dG3D/benchmarks/interpolationPBC_2DShifted/idealHole.py b/dG3D/benchmarks/interpolationPBC_2DShifted/idealHole.py
index cccb02991f039b7fdfc00be342b3c2405e13d3a2..74b3bd15bd4adef6b580e93aace20a78253a7fe9 100644
--- a/dG3D/benchmarks/interpolationPBC_2DShifted/idealHole.py
+++ b/dG3D/benchmarks/interpolationPBC_2DShifted/idealHole.py
@@ -64,17 +64,15 @@ mysolver.setPeriodicity(0,1.,0,"y")
mysolver.setPeriodicity(0,0,1.,"z")
#boundary condition
-microBC = nonLinearPeriodicBC(100,2)
+microBC = nonLinearShiftedPeriodicBC(100,2)
microBC.setOrder(1)
microBC.setBCPhysical(1,4,3,2)
method =1 # Periodic mesh = 0, Langrange interpolation = 1, Cubic spline interpolation =2, FE linear= 3, FE Quad = 4
degree = 5 # Order used for polynomial interpolation
-addvertex = 0 # Polynomial interpolation by mesh vertex = 0, Polynomial interpolation by virtual vertex
+addvertex = 1 # Polynomial interpolation by mesh vertex = 0, Polynomial interpolation by virtual vertex
microBC.setPeriodicBCOptions(method, degree,bool(addvertex))
theta = pi/10;
-
-microBC.setShiftedBC(True)
microBC.setShiftPBCNormal(sin(theta),-cos(theta),0.)
# Deformation gradient
diff --git a/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_changeBCAfterFailure/model.py b/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_changeBCAfterFailure/model.py
index 5a8cdc246c94f3bda8e931cfe709c7037f49c68f..baee0521ae14f52f9b72bedc7b06c9a96845d2f1 100644
--- a/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_changeBCAfterFailure/model.py
+++ b/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_changeBCAfterFailure/model.py
@@ -100,7 +100,7 @@ lawCoh = dG3DMultiscaleCohesiveLaw(lcohNum,0.7)
lawCoh.setCharacteristicLength(1.)
lawCoh.setLostSolutionUniquenssTolerance(0.1)
-microBCFailure = nonLinearPeriodicBC(1000,2)
+microBCFailure = nonLinearShiftedPeriodicBC(1000,2)
microBCFailure.setOrder(1)
microBCFailure.setBCPhysical(1,4,3,2)
# periodiodic BC
@@ -108,7 +108,6 @@ method = 1 # Periodic mesh = 0 Langrange interpolation = 1 Cubic spline interpol
degree = 5
addvertex = False
microBCFailure.setPeriodicBCOptions(method, degree,addvertex)
-microBCFailure.setShiftedBC(True) # shifted BC at interface only
lawCoh.addMicroBCFailure(microBCFailure)
diff --git a/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_shiftedBC/model.py b/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_shiftedBC/model.py
index 88b0c000f402c1c1fc0d1d6b1adc112949a99c14..5cda30b84f825ca75c8c9d127f9b7abceaeda837 100644
--- a/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_shiftedBC/model.py
+++ b/dG3D/benchmarks/multiscaleCohesiveTest2D_fullDG_shiftedBC/model.py
@@ -29,7 +29,7 @@ micromeshfile="micro.msh" # name of mesh file
nfield = 11 # number of the field (physical number of entity)
myfield1 = nonLocalDamageDG3DDomain(1000,nfield,0,lawnum1,0,1e3,2,1)
-microBC = nonLinearPeriodicBC(1000,2)
+microBC = nonLinearShiftedPeriodicBC(1000,2)
microBC.setOrder(1)
microBC.setBCPhysical(1,4,3,2)
# periodiodic BC
@@ -37,7 +37,6 @@ method = 1 # Periodic mesh = 0 Langrange interpolation = 1 Cubic spline interpol
degree = 5
addvertex = False
microBC.setPeriodicBCOptions(method, degree,addvertex)
-microBC.setShiftedBC(True) # shifted BC at interface only
# DEFINE MACROPROBLEM
matnum1 = 1;
diff --git a/dG3D/src/dG3DMultiscaleMaterialLaw.cpp b/dG3D/src/dG3DMultiscaleMaterialLaw.cpp
index fe8c29e000f91a39ba4d9629f9f829ac93116395..24c1231a5efc4b3df2d5a2916fdfb70e3bf6df59 100644
--- a/dG3D/src/dG3DMultiscaleMaterialLaw.cpp
+++ b/dG3D/src/dG3DMultiscaleMaterialLaw.cpp
@@ -380,12 +380,10 @@ void MultiscaleFractureByCohesive3DLaw::createIPState(const bool isSolve, IPStat
mixedBC->setKUBCDirection(b);
}
}
- else if (mbc->getType() == nonLinearMicroBC::PBC){
- nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(mbc);
+ else if (mbc->getType() == nonLinearMicroBC::ShiftedPBC){
+ nonLinearShiftedPeriodicBC* pbc = static_cast<nonLinearShiftedPeriodicBC*>(mbc);
if (mbc->getDim() == 2){
- if (pbc->isShifted()){
- pbc->setShiftPBCNormal(n);
- }
+ pbc->setShiftPBCNormal(n);
}
else{
Msg::Fatal("shifted BC is not defined for 3D problem");