diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
index 3dfe11b0f7f45be084c5a0ac3c7322b4296b8d7c..4e1e963974ff1024cc565000a31b80b7c90b075d 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.cpp
@@ -278,18 +278,9 @@ void nonLinearPeriodicBC::setPeriodicBCOptions(const int method, const int seg,
this->setPBCPolynomialDegree(1,seg);
this->setPBCPolynomialDegree(2,seg);
- if (method == 3 or method == 4){
- printf("Using FE segment method, add new vertex flag is always actived \n");
- this->setAddNewVertices(0,true);
- this->setAddNewVertices(1,true);
- this->setAddNewVertices(2,true);
- }
- else{
- this->setAddNewVertices(0,add);
- this->setAddNewVertices(1,add);
- this->setAddNewVertices(2,add);
- }
-
+ this->setAddNewVertices(0,add);
+ this->setAddNewVertices(1,add);
+ this->setAddNewVertices(2,add);
if (this->getPBCMethod() == nonLinearPeriodicBC::CEM){
printf("Periodic mesh formulation is used for building the PBC constraint \n");
diff --git a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
index 52dab7f328d3a408a98469c1265c2a4d06d26358..c7cfcc42481405439862b5191eabdf548b669d84 100644
--- a/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
+++ b/NonLinearSolver/BoundaryConditions/nonLinearMicroBC.h
@@ -442,18 +442,11 @@ class nonLinearPeriodicBC: public nonLinearMicroBC{
return _forceInterp;
}
- bool isNewVertices() const{
- return (_XaddNewVertices&& _YaddNewVertices)&&_ZaddNewVertices;
- }
-
bool addNewVertices(const int i) const { // add new Dof flag
- if (_wM == CSIM) return true;
- else {
- if (i ==0) return _XaddNewVertices;
- else if (i==1) return _YaddNewVertices;
- else if (i==2) return _ZaddNewVertices;
- else Msg::Fatal("direction must be correctly defined addNewVertices(i)");
- };
+ if (i ==0) return _XaddNewVertices;
+ else if (i==1) return _YaddNewVertices;
+ else if (i==2) return _ZaddNewVertices;
+ else Msg::Fatal("direction must be correctly defined addNewVertices(i)");
};
whichMethod getPBCMethod() const{
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index dc33c05dac037bc13832510b3cd0ae6b7c52693c..7a3cc2d27db5e7179a4e9611e5455a67b2dde87e 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -9051,8 +9051,7 @@ stiffnessCondensation* nonLinearMechSolver::createStiffnessCondensation(const bo
}
else if (dynamic_cast<nonLinearPeriodicBC*>(_microBC)){
nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_microBC);
- bool addnewvertexflag = pbc->isNewVertices();
- if (addnewvertexflag and pbc->getOrder()==1){
+ if (_pAl->getSplittedDof()->sizeOfVirtualDof() >0 and pbc->getOrder()==1){
#if defined(HAVE_PETSC)
Msg::Info("interpolationStiffnessCondensationPETSC is used for nonLinearPeriodicBC");
stiffCon = new interpolationStiffnessCondensationPETSC(pAssembler,_pAl,sflag,tflag);
diff --git a/NonLinearSolver/periodicBC/numericalFunctions.h b/NonLinearSolver/periodicBC/numericalFunctions.h
index 8908cbf961e6cd65cfd1482a80acd033eb82224a..559a2eaf4c299823baf39ca6d281542d7d6f5cb5 100644
--- a/NonLinearSolver/periodicBC/numericalFunctions.h
+++ b/NonLinearSolver/periodicBC/numericalFunctions.h
@@ -19,6 +19,22 @@
#include "STensor53.h"
#include "STensor63.h"
+inline SPoint3 computeTwoLineIntersectionOnSurface(const SPoint3& A, const SPoint3& B, const SPoint3& C, const SPoint3& D){
+ // Line AB and CD cut at I
+ // ABCD on a same plane
+ SVector3 AB(A,B);
+ SVector3 CD(C,D);
+ SVector3 facN = crossprod(AB,CD);
+ SVector3 n = crossprod(CD,facN);
+ SVector3 AD(A,D);
+ double t = dot(AD,n)/dot(AB,n);
+
+ double x = A.x() + t*AB.x();
+ double y = A.y() + t*AB.y();
+ double z = A.z() + t*AB.z();
+ return SPoint3(x,y,z);
+};
+
inline void computeSurface(const SPoint3& A, const SPoint3& B, const SPoint3& C,double &a, double& b, double &c, double &d){
SVector3 AB(A,B);
SVector3 AC(A,C);
diff --git a/NonLinearSolver/periodicBC/pbcAlgorithm.cpp b/NonLinearSolver/periodicBC/pbcAlgorithm.cpp
index c8e943641e10154aec7fd39aa3dbef1d7b02fbf6..2fb6e3b6d853eafbe4df44e6c36ffe9883761d56 100644
--- a/NonLinearSolver/periodicBC/pbcAlgorithm.cpp
+++ b/NonLinearSolver/periodicBC/pbcAlgorithm.cpp
@@ -120,13 +120,6 @@ int pbcAlgorithm::getNumberMicroToMacroHomogenizedVariables() const{
return dimStress;
}
-void pbcAlgorithm::switchMicroBC(){
- _pbcConstraint->switchMicroBC();
- // split dofs
- this->splitDofs();
-};
-
-
void pbcAlgorithm::updateConstraint(const IPField* ipfield){
for (std::set<constraintElement*>::iterator it= _pbcConstraint->constraintGroupBegin();
it!= _pbcConstraint->constraintGroupEnd(); it++){
diff --git a/NonLinearSolver/periodicBC/pbcAlgorithm.h b/NonLinearSolver/periodicBC/pbcAlgorithm.h
index 5e56af6eb1fe0110085854babdb5def04521ad0f..e9ed9519f781d82224e76db36c73b0dbd2d1dba6 100644
--- a/NonLinearSolver/periodicBC/pbcAlgorithm.h
+++ b/NonLinearSolver/periodicBC/pbcAlgorithm.h
@@ -36,8 +36,6 @@ class pbcAlgorithm{
pbcAlgorithm(nonLinearMechSolver* solver);
virtual ~pbcAlgorithm();
- void switchMicroBC();
-
pbcDofManager* getSplittedDof() {return _splittedDofs;};
splitStiffnessMatrix* getSplitStiffness() {return _splittedDofs->getSplitStiffnessMatrix();};
pbcConstraintElementGroup* getPBCConstraintGroup() {return _pbcConstraint;};
diff --git a/NonLinearSolver/periodicBC/pbcConstraintElement.cpp b/NonLinearSolver/periodicBC/pbcConstraintElement.cpp
index 129a74395f77914abc4bf878161c98ef398b8ea8..f84b5cc020cd10fe983aad60245142654730f8ec 100644
--- a/NonLinearSolver/periodicBC/pbcConstraintElement.cpp
+++ b/NonLinearSolver/periodicBC/pbcConstraintElement.cpp
@@ -631,8 +631,9 @@ FEConstraintElement::FEConstraintElement(nonLinearMicroBC* mbc, FunctionSpaceBas
};
_Cbc.resize(_ndofs,size*_ndofs);
_Cbc.setAll(0.);
- for (int j=0; j<_ndofs; j++){
- for (int i=0; i<size; i++){
+
+ for (int i=0; i<size; i++){
+ for (int j=0; j<_ndofs; j++){
_C(j,col) = -1.*_fVal[i];
_Cbc(j,col-_ndofs) = _fVal[i];
col++;
@@ -672,7 +673,9 @@ FEConstraintElement::~FEConstraintElement(){};
void FEConstraintElement::getConstraintKeys(std::vector<Dof>& keys) const{
getKeysFromVertex(_periodicSpace,_v, getComp(), keys);
- _periodicSpace->getKeys(_e,keys);
+ for (int i=0; i< _e->getNumVertices(); i++){
+ getKeysFromVertex(_periodicSpace,_e->getVertex(i),getComp(),keys);
+ }
};
void FEConstraintElement::getMultiplierKeys(std::vector<Dof>& keys) const{
diff --git a/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp b/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
index 9b9bb338cb0f928de52000f35c7d39d596526ec2..816d95a90707485c8eeee8e00d3d6411f47b1bb1 100644
--- a/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
+++ b/NonLinearSolver/periodicBC/pbcCreateConstraints.cpp
@@ -19,6 +19,374 @@
#include "pbcSupplementConstraintExtraDofDiffusion.h"
#include "directionalConstraintElement.h"
+void InterpolationOperations::groupIntersection(groupOfElements* g1, groupOfElements* g2, std::set<MVertex*>& ver){
+ for (groupOfElements::vertexContainer::iterator it1 = g1->vbegin(); it1!= g1->vend(); it1++){
+ MVertex* v1 = *it1;
+ for (groupOfElements::vertexContainer::iterator it2 = g2->vbegin(); it2!= g2->vend(); it2++){
+ MVertex* v2 = *it2;
+ if (v1->getNum() ==v2->getNum()){
+ ver.insert(v1);
+ };
+ };
+ };
+};
+
+void InterpolationOperations::groupIntersection(std::set<MVertex*>& g1, std::set<MVertex*>& g2, std::set<MVertex*>& ver){
+ for (std::set<MVertex*>::iterator it1 = g1.begin(); it1!= g1.end(); it1++){
+ MVertex* v1 = *it1;
+ for (std::set<MVertex*>::iterator it2 = g2.begin(); it2 != g2.end(); it2++){
+ MVertex* v2 = *it2;
+ if (v1==v2){
+ ver.insert(v1);
+ };
+ };
+ };
+};
+
+void InterpolationOperations::intersectionBetweenTwoFaces(const groupOfElements* face1, const groupOfElements* face2, groupOfElements& line){
+
+ std::vector<std::vector<MVertex*> > allEdges;
+ for (groupOfElements::elementContainer::iterator itL = face1->begin(); itL!= face1->end(); itL++){
+ MElement* eleL = *itL;
+ if (eleL->getDim() != 2){
+ Msg::Fatal("This fuction is implemented from 2D elements only");
+ }
+ bool found = false;
+ for (int i=0; i< eleL->getNumEdges(); i++){
+ std::vector<MVertex*> vvL;
+ eleL->getEdgeVertices(i,vvL);
+ // found in face3
+ for (groupOfElements::elementContainer::iterator itR = face2->begin(); itR!= face2->end(); itR++){
+ MElement* eleR = *itR;
+ for (int j=0; j< eleR->getNumEdges(); j++){
+ std::vector<MVertex*> vvR;
+ eleR->getEdgeVertices(j,vvR);
+ if (((vvL[0]->getNum() == vvR[0]->getNum()) and (vvL[1]->getNum() == vvR[1]->getNum())) or
+ ((vvL[0]->getNum() == vvR[1]->getNum()) and (vvL[1]->getNum() == vvR[0]->getNum()))){
+ allEdges.push_back(vvL);
+ found = true;
+ break;
+ }
+
+ }
+ if (found) break;
+ }
+ if (found) break;
+
+ }
+ }
+
+ // add to grEle
+ for (int i=0; i< allEdges.size(); i++){
+ std::vector<MVertex*>& vv = allEdges[i];
+ MElement* ele = new MLineN(vv);
+ line.insert(ele);
+ }
+};
+
+void InterpolationOperations::createLinePolynomialBase(const int degree, // degree
+ const bool createNewVer, // true if new vertices are created false if interpolation base is taken from existing vertices
+ MVertex* vstart, MVertex* vend, // vstart and vend
+ std::set<MVertex*>& line, // current line
+ std::vector<MVertex*>& base, // base obtained
+ std::set<MVertex*>& virtVer // virtualVertices to store all virtual virtices which has been created
+ ){
+ base.resize(degree+1);
+ base[0] = vstart;
+ base[degree] = vend;
+ if (degree > 1){
+ if (createNewVer){
+ // create (degree-1) new vertices which are equally distributed
+ double x0 = vstart->x(); double x1 = vend->x();
+ double y0 = vstart->y(); double y1 = vend->y();
+ double z0 = vstart->z(); double z1 = vend->z();
+
+ double dx = (x1-x0)/degree;
+ double dy = (y1-y0)/degree;
+ double dz = (z1-z0)/degree;
+ for (int j=0; j<degree-1; j++){
+ double xi = x0+ (j+1)*dx;
+ double yi = y0+ (j+1)*dy;
+ double zi = z0+ (j+1)*dz;
+ MVertex* v = new MVertex(xi,yi,zi);
+ base[j+1] = v;
+ virtVer.insert(v); // as new vertices
+ };
+ }
+ else{
+ // get (degree-1) vertices from existing vertices
+ std::vector<MVertex*> list;
+ std::vector<double> vals;
+ for (std::set<MVertex*>::iterator it=line.begin(); it!=line.end(); it++){
+ MVertex* v=*it;
+ if ((v->getNum() != vstart->getNum()) and (v->getNum() != vend->getNum())){
+ list.push_back(v);
+ vals.push_back(v->distance(vstart));
+ }
+ };
+
+ // reorder based on distant to vstart
+ int size=list.size();
+ for (int i=0; i<size; i++){
+ for (int j=i; j<size; j++){
+ if (vals[j]<vals[i]){
+ double temp=vals[i]; vals[i]=vals[j]; vals[j]=temp;
+ MVertex* vtemp=list[i]; list[i]=list[j]; list[j]=vtemp;
+ };
+ };
+ };
+ // at leat
+ if (degree-1>=size){
+ Msg::Warning("all vertices on this linea are taken in interpolation basis");
+ for (int i=0; i< size; i++){
+ base[i+1] = list[i];
+ }
+ }
+ else{
+ int interval = (size)/(degree-1);
+ for (int j=0; j<degree-1; j++){
+ base[j+1] = list[j*interval];
+ };
+ }
+ }
+ }
+};
+
+
+template<class Iterator>
+MVertex* InterpolationOperations::findNearestVertex(Iterator itbegin, Iterator itend, MVertex* vp){
+ MVertex* vinit = *itbegin;
+ double dist=vinit->distance(vp);
+ for (Iterator it=itbegin; it!=itend; it++){
+ MVertex* v=*it;
+ double tempDist=v->distance(vp);
+ if (tempDist<dist) {
+ dist=tempDist;
+ vinit=v;
+ };
+ };
+ return vinit;
+};
+
+MVertex* InterpolationOperations::findNearestVertex(groupOfElements* g, MVertex* vp){
+
+ return findNearestVertex(g->vbegin(), g->vend(),vp);
+};
+
+void InterpolationOperations::createGroupLinearLineElementFromVerticesVector(std::vector<MVertex*>& b, groupOfElements& gr){
+ for (int i=0; i<b.size()-1; i++){
+ MElement* line = new MLine(b[i],b[i+1]);
+ gr.insert(line);
+ }
+};
+
+void InterpolationOperations::createGroupQuadraticLineElementFromVerticesVector(std::vector<MVertex*>& b, groupOfElements& gr, std::set<MVertex*>& virVertices){
+ // base will be modified due to create new vertices
+ std::vector<MVertex*> xtemp(b);
+ b.clear();
+ for (int i=0; i<xtemp.size()-1; i++){
+ b.push_back(xtemp[i]);
+ MVertex* v = new MVertex(0.5*(xtemp[i+1]->x()+xtemp[i]->x()),
+ 0.5*(xtemp[i+1]->y()+xtemp[i]->y()),
+ 0.5*(xtemp[i+1]->z()+xtemp[i]->z()));
+
+ b.push_back(v);
+ virVertices.insert(v);
+ MElement* line = new MLine3(xtemp[i],xtemp[i+1],v);
+ gr.insert(line);
+ }
+ // add last vertex
+ b.push_back(xtemp[xtemp.size()-1]);
+};
+
+void InterpolationOperations::createLinearElementBCGroup3DFromTwoEdges( std::vector<MVertex*>& vlx, std::vector<MVertex*>& vly,
+ MVertex* vc, groupOfElements& g, std::set<MVertex*>& newVer){
+ if (vlx[0] != vly[0]){
+ Msg::Fatal("two BC group must be started by same vertex");
+ return;
+ }
+ MVertex* vroot = vlx[0];
+
+ int NX = vlx.size();
+ int NY = vly.size();
+
+ std::vector<std::vector<MVertex*> > vmat;
+ vmat.push_back(vlx);
+
+ for (int j=1; j<NY; j++){
+ std::vector<MVertex*> vv;
+ vv.push_back(vly[j]);
+ for (int i=1; i< NX; i++){
+ MVertex* vertex = NULL;
+ if (j == NY-1 && i == NX-1){
+ vertex = vc;
+ }
+ else{
+ vertex = new MVertex(0,0,0);
+ vertex->x() = vlx[i]->x() - vroot->x()+ vly[j]->x();
+ vertex->y() = vlx[i]->y() - vroot->y()+ vly[j]->y();
+ vertex->z() = vlx[i]->z() - vroot->z()+ vly[j]->z();
+ newVer.insert(vertex);
+ }
+ vv.push_back(vertex);
+ }
+ vmat.push_back(vv);
+ }
+
+
+ for (int i=0; i<NX-1; i++){
+ for (int j=0; j<NY-1; j++){
+ MElement* e = new MQuadrangle(vmat[j][i],vmat[j][i+1],vmat[j+1][i+1],vmat[j+1][i]);
+ g.insert(e);
+ }
+ }
+
+};
+
+void InterpolationOperations::createQuadraticElementBCGroup3DFromTwoEdges( std::vector<MVertex*>& vlx, std::vector<MVertex*>& vly,
+ MVertex* vc, groupOfElements& g, std::set<MVertex*>& newVer){
+ if (vlx[0] != vly[0]){
+ Msg::Fatal("two BC group must be started by same vertex");
+ return;
+ }
+ MVertex* vroot = vlx[0];
+
+ int NX = vlx.size();
+ int NY = vly.size();
+
+ std::vector<std::vector<MVertex*> > vmat;
+ vmat.push_back(vlx);
+
+ for (int j=1; j<NY; j++){
+ std::vector<MVertex*> vv;
+ vv.push_back(vly[j]);
+ for (int i=1; i< NX; i++){
+ MVertex* vertex = NULL;
+ if (j == NY-1 && i == NX-1){
+ vertex = vc;
+ }
+ else{
+ vertex = new MVertex(0,0,0);
+ vertex->x() = vlx[i]->x() - vroot->x()+ vly[j]->x();
+ vertex->y() = vlx[i]->y() - vroot->y()+ vly[j]->y();
+ vertex->z() = vlx[i]->z() - vroot->z()+ vly[j]->z();
+ newVer.insert(vertex);
+ }
+ vv.push_back(vertex);
+ }
+ vmat.push_back(vv);
+ }
+
+
+ for (int i=0; i<NX-1; i+=2){
+ for (int j=0; j<NY-1; j+=2){
+ MElement* e = new MQuadrangle9(vmat[j][i],vmat[j][i+2],vmat[j+2][i+2],vmat[j+2][i],
+ vmat[j][i+1],vmat[j+1][i+2],vmat[j+2][i+1],vmat[j+1][i],
+ vmat[j+1][i+1]);
+ g.insert(e);
+ }
+ }
+
+};
+
+void InterpolationOperations::createPolynomialBasePoints(nonLinearPeriodicBC* pbc,
+ MVertex* v1, MVertex* v2, MVertex* v4, MVertex* v5,
+ std::set<MVertex*>& l12, std::set<MVertex*>& l41,std::set<MVertex*>& l15,
+ std::vector<MVertex*>& xbase, std::vector<MVertex*>& ybase, std::vector<MVertex*>& zbase,
+ std::set<MVertex*>& newVertices){
+ // direction x, y, z are defined as follows
+ // x - 12
+ // y - 14
+ // z - 15
+
+ int dim = pbc->getDim();
+
+ int degreeX = pbc->getPBCPolynomialDegree(0);
+ int degreeY = pbc->getPBCPolynomialDegree(1);;
+ int degreeZ = pbc->getPBCPolynomialDegree(2);;
+
+ // modify degree following available vertices in l12, l14, l15
+ if (!pbc->isForcedInterpolationDegree()){
+ int sizeX = l12.size();
+ int sizeY = l41.size();
+ int sizeZ = l15.size();
+ Msg::Info("Maximal number of Boundary Nodes in X direction %d, Y direction %d, Z direction %d",sizeX,sizeY,sizeZ);
+
+
+ if (pbc->getPBCMethod() == nonLinearPeriodicBC::LIM or
+ pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN){
+ // first order
+ if (degreeX>sizeX-1) degreeX = sizeX-1;
+ if (degreeY>sizeY-1) degreeY = sizeY-1;
+ if (degreeZ>sizeZ-1 && dim == 3) degreeZ = sizeZ-1;
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM or
+ pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
+ // second-order
+ if (degreeX>(sizeX-1)/2) degreeX = (sizeX-1)/2;
+ if (degreeY>(sizeY-1)/2) degreeY = (sizeY-1)/2;
+ if (degreeZ>(sizeZ-1)/2 && dim == 3) degreeZ = (sizeZ-1)/2;
+ }
+ else{
+ Msg::Fatal("interpolation type must be defined");
+ }
+ // at leat, interpolation must be linear
+ if (degreeX<1) degreeX = 1;
+ if (degreeY<1) degreeY = 1;
+ if (degreeZ<1) degreeZ = 1;
+
+ // modify degree
+ pbc->setPBCPolynomialDegree(0,degreeX);
+ pbc->setPBCPolynomialDegree(1,degreeY);
+ pbc->setPBCPolynomialDegree(2,degreeZ);
+
+ // force add new vertices
+ if (sizeX < 2) pbc->setAddNewVertices(0,true);
+ if (sizeY < 2) pbc->setAddNewVertices(1,true);
+ if (sizeZ < 2) pbc->setAddNewVertices(2,true);
+
+ }
+
+ Msg::Info("Interpolation degree used in X direction %d, Y direction %d, Z direction %d",degreeX,degreeY,degreeZ);
+
+
+ InterpolationOperations::createLinePolynomialBase(degreeX,pbc->addNewVertices(0),v1,v2,l12,xbase,newVertices);
+ InterpolationOperations::createLinePolynomialBase(degreeY,pbc->addNewVertices(1),v1,v4,l41,ybase,newVertices);
+ if (dim == 3){
+ InterpolationOperations::createLinePolynomialBase(degreeZ,pbc->addNewVertices(2),v1,v5,l15,zbase,newVertices);
+ }
+
+
+};
+
+bool InterpolationOperations::isInside(MVertex* v, MElement* e){
+ int dim = e->getDim();
+ std::vector<MVertex*> vv;
+ e->getVertices(vv);
+
+ SPoint3 pp;
+ if (dim == 1){
+ pp = project(v->point(),vv[0]->point(), vv[1]->point());
+ }
+ else if (dim ==2){
+ pp = project(v->point(),vv[0]->point(), vv[1]->point(),vv[2]->point());
+ }
+
+ double xyz[3] = {pp[0],pp[1],pp[2]};
+ double uvw[3];
+
+ e->xyz2uvw(xyz,uvw);
+ if (e->isInside(uvw[0],uvw[1],uvw[2])){
+ return true;
+ }
+ else return false;
+
+};
+
+
+
+
+
groupOfElementsDomainDecomposition::groupOfElementsDomainDecomposition(const int phys, groupOfElements* g_, std::vector<partDomain*>& allDomain):
g(g_),physical(phys){
// init map
@@ -77,12 +445,12 @@ groupOfElementsDomainDecomposition::~groupOfElementsDomainDecomposition(){
pbcConstraintElementGroup::pbcConstraintElementGroup(nonLinearMechSolver* sl,
FunctionSpaceBase* lagspace, FunctionSpaceBase* multspace)
:_solver(sl), _v1(NULL),_v2(NULL),_v4(NULL),_v5(NULL),
- _v3(NULL),_v6(NULL),_v7(NULL),_v8(NULL), _gXY(NULL),_gYZ(NULL),_gZX(NULL),
- _gX(NULL),_gY(NULL),_gZ(NULL), _g(NULL),
+ _v3(NULL),_v6(NULL),_v7(NULL),_v8(NULL), _g(NULL),
_LagSpace(lagspace),_MultSpace(multspace),
_periodic12(0.),_periodic14(0.),_periodic15(0.),
- _e1(1.,0.,0.),_e2(0.,1.,0.),_e3(0.,0.,1.),
- _e1Rotated(1.,0.,0.),_e2Rotated(0.,1.,0.),_e3Rotated(0.,0.,1.){
+ f1234(NULL), f5678(NULL), f1562(NULL),f4873(NULL), f1584(NULL), f2673(NULL),
+ _gX(NULL),_gY(NULL),_gZ(NULL),
+ _gXY(NULL),_gYZ(NULL),_gZX(NULL){
this->__init__();
};
@@ -92,82 +460,29 @@ pbcConstraintElementGroup::~pbcConstraintElementGroup(){
if (cel) delete cel;
cel = NULL;
}
- _allConstraint.clear();
- _cVertices.clear();
- _bVertices.clear();
for (std::set<MVertex*>::iterator it = _virtualVertices.begin(); it!= _virtualVertices.end(); it++){
MVertex* v = *it;
if (v) delete v;
v = NULL;
}
- _virtualVertices.clear();
- l12.clear();
- l23.clear();
- l34.clear();
- l41.clear();
- l56.clear();
- l67.clear();
- l78.clear();
- l85.clear();
- l15.clear();
- l26.clear();
- l37.clear();
- l48.clear();
if (_g) delete _g; _g = NULL;
- if (_gX) delete _gX; _gX = NULL;
- if (_gY) delete _gY; _gY = NULL;
- if (_gZ) delete _gZ; _gZ = NULL;
-
- if (_gXY) delete _gXY; _gXY = NULL;
- if (_gYZ) delete _gYZ; _gYZ = NULL;
- if (_gZX) delete _gZX; _gZX = NULL;
+
+ if (_gX) delete _gX;
+ if (_gY) delete _gY;
+ if (_gZ) delete _gZ;
+ if (_gXY) delete _gXY;
+ if (_gYZ) delete _gYZ;
+ if (_gZX) delete _gZX;
for (int i=0; i< _boundaryDGMap.size(); i++){
delete _boundaryDGMap[i];
}
- _boundaryDGMap.clear();
for (std::map<partDomain*,FunctionSpaceBase*>::iterator it = _allDGMultSpace.begin(); it!=_allDGMultSpace.end(); it++){
if (it->second!=NULL) delete it->second;
}
- _allDGMultSpace.clear();
};
-void pbcConstraintElementGroup::groupIntersection(groupOfElements* g1, groupOfElements* g2, std::set<MVertex*>& ver) const{
- for (groupOfElements::vertexContainer::iterator it1 = g1->vbegin(); it1!= g1->vend(); it1++){
- MVertex* v1 = *it1;
- for (groupOfElements::vertexContainer::iterator it2 = g2->vbegin(); it2!= g2->vend(); it2++){
- MVertex* v2 = *it2;
- if (v1->getNum() ==v2->getNum()){
- ver.insert(v1);
- };
- };
- };
-};
-void pbcConstraintElementGroup::groupIntersection(groupOfElements* g1, groupOfElements* g2, groupOfElements* &g) const{
- std::vector<MElement*> egroup;
- for (groupOfElements::elementContainer::iterator it1 = g1->begin(); it1!= g1->end(); it1++){
- MElement* e1 = *it1;
- for (groupOfElements::elementContainer::iterator it2 = g2->begin(); it2!= g2->end(); it2++){
- MElement* e2 = *it2;
- if (e1 == e2){
- egroup.push_back(e1);
- };
- };
- };
- g = new groupOfElements(egroup);
-};
-void pbcConstraintElementGroup::groupIntersection(std::set<MVertex*>& g1, std::set<MVertex*>& g2, std::set<MVertex*>& ver) const{
- for (std::set<MVertex*>::iterator it1 = g1.begin(); it1!= g1.end(); it1++){
- MVertex* v1 = *it1;
- for (std::set<MVertex*>::iterator it2 = g2.begin(); it2 != g2.end(); it2++){
- MVertex* v2 = *it2;
- if (v1==v2){
- ver.insert(v1);
- };
- };
- };
-};
void pbcConstraintElementGroup::__init__(){
std::vector<partDomain*>& allDomain = *(_solver->getDomainVector());
@@ -193,12 +508,16 @@ void pbcConstraintElementGroup::__init__(){
if (_g->size()==0)
Msg::Fatal("Error in setting domain");
+
+ const std::vector<int>& bphysical = _solver->getMicroBC()->getBoundaryPhysicals();
+ const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
+ if (bgroup.size() == 0) Msg::Fatal("microscopic boundary condition has not any support");
+ int dim = _solver->getMicroBC()->getDim();
if (_fullDG){
Msg::Info("start decomsing boundary element to domains");
_boundaryDGMap.clear();
- const std::vector<int>& bphysical = _solver->getMicroBC()->getBoundaryPhysicals();
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
+
for (int ig=0; ig< bgroup.size(); ig++){
groupOfElements* gr = bgroup[ig];
groupOfElementsDomainDecomposition* grDomPart = new groupOfElementsDomainDecomposition(bphysical[ig], gr,allDomain);
@@ -227,361 +546,235 @@ void pbcConstraintElementGroup::__init__(){
}
}
}
+
+ if (dim == 2){
+ if (bgroup.size() != 4) Msg::Fatal("Boundary support must be 4 edges");
+ l12.clear();
+ l12.insert(bgroup[0]->vbegin(),bgroup[0]->vend());
+ l41.clear();
+ l41.insert(bgroup[1]->vbegin(),bgroup[1]->vend());
+ l34.clear();
+ l34.insert(bgroup[2]->vbegin(),bgroup[2]->vend());
+ l23.clear();
+ l23.insert(bgroup[3]->vbegin(),bgroup[3]->vend());
+
+ _bVertices.clear();
+ InterpolationOperations::groupIntersection(bgroup[0],bgroup[1],_bVertices);
+ InterpolationOperations::groupIntersection(bgroup[1],bgroup[2],_bVertices);
+ InterpolationOperations::groupIntersection(bgroup[2],bgroup[3],_bVertices);
+ InterpolationOperations::groupIntersection(bgroup[3],bgroup[0],_bVertices);
+
+ _cVertices.clear();
+ _cVertices.insert(_bVertices.begin(),_bVertices.end());
+ #ifdef _DEBUG
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ Msg::Info("corner vertex %d ",(*it)->getNum());
+ };
+ #endif //_DEBUG
+ }
+ else if (dim == 3){
+ if (bgroup.size() != 6) Msg::Fatal("Boundary support must be 6 faces");
+ // all faces
+ f1234 = bgroup[0];
+ f5678 = bgroup[3];
+ f1562 = bgroup[1];
+ f4873 = bgroup[4];
+ f1584 = bgroup[2];
+ f2673 = bgroup[5];
+
+ // all edges
+ l12.clear(); InterpolationOperations::groupIntersection(f1234,f1562,l12);
+ l23.clear(); InterpolationOperations::groupIntersection(f1234,f2673,l23);
+ l34.clear(); InterpolationOperations::groupIntersection(f1234,f4873,l34);
+ l41.clear(); InterpolationOperations::groupIntersection(f1234,f1584,l41);
+ l56.clear(); InterpolationOperations::groupIntersection(f5678,f1562,l56);
+ l67.clear(); InterpolationOperations::groupIntersection(f5678,f2673,l67);
+ l78.clear(); InterpolationOperations::groupIntersection(f5678,f4873,l78);
+ l85.clear(); InterpolationOperations::groupIntersection(f5678,f1584,l85);
+ l15.clear(); InterpolationOperations::groupIntersection(f1562,f1584,l15);
+ l26.clear(); InterpolationOperations::groupIntersection(f1562,f2673,l26);
+ l37.clear(); InterpolationOperations::groupIntersection(f2673,f4873,l37);
+ l48.clear(); InterpolationOperations::groupIntersection(f4873,f1584,l48);
+
+
+ _bVertices.clear();
+ _bVertices.insert(l12.begin(),l12.end());
+ _bVertices.insert(l23.begin(),l23.end());
+ _bVertices.insert(l34.begin(),l34.end());
+ _bVertices.insert(l41.begin(),l41.end());
+ _bVertices.insert(l56.begin(),l56.end());
+ _bVertices.insert(l67.begin(),l67.end());
+ _bVertices.insert(l78.begin(),l78.end());
+ _bVertices.insert(l85.begin(),l85.end());
+ _bVertices.insert(l15.begin(),l15.end());
+ _bVertices.insert(l26.begin(),l26.end());
+ _bVertices.insert(l37.begin(),l37.end());
+ _bVertices.insert(l48.begin(),l48.end());
+
+
+ #ifdef _DEBUG
+ for (std::set<MVertex*>::iterator it = l12.begin(); it!= l12.end(); it++){
+ Msg::Info("vertex %d 12",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l23.begin(); it!= l23.end(); it++){
+ Msg::Info("vertex %d 23",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l34.begin(); it!= l34.end(); it++){
+ Msg::Info("vertex %d 34",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l41.begin(); it!= l41.end(); it++){
+ Msg::Info("vertex %d 41",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l56.begin(); it!= l56.end(); it++){
+ Msg::Info("vertex %d 56",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l67.begin(); it!= l67.end(); it++){
+ Msg::Info("vertex %d 67",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l78.begin(); it!= l78.end(); it++){
+ Msg::Info("vertex %d 78",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l85.begin(); it!= l85.end(); it++){
+ Msg::Info("vertex %d 85",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l15.begin(); it!= l15.end(); it++){
+ Msg::Info("vertex %d 15",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l26.begin(); it!= l26.end(); it++){
+ Msg::Info("vertex %d 26",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l37.begin(); it!= l37.end(); it++){
+ Msg::Info("vertex %d 37",(*it)->getNum());
+ };
+ for (std::set<MVertex*>::iterator it = l48.begin(); it!= l48.end(); it++){
+ Msg::Info("vertex %d 48",(*it)->getNum());
+ };
+ #endif //_DEBUG
+
+ // corner vertices
+ _cVertices.clear();
+ InterpolationOperations::groupIntersection(l12,l23,_cVertices);
+ InterpolationOperations::groupIntersection(l23,l34,_cVertices);
+ InterpolationOperations::groupIntersection(l34,l41,_cVertices);
+ InterpolationOperations::groupIntersection(l41,l12,_cVertices);
+ InterpolationOperations::groupIntersection(l56,l67,_cVertices);
+ InterpolationOperations::groupIntersection(l67,l78,_cVertices);
+ InterpolationOperations::groupIntersection(l78,l85,_cVertices);
+ InterpolationOperations::groupIntersection(l85,l56,_cVertices);
+
+ #ifdef _DEBUG
+ for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ Msg::Info("corner vertex %d ",(*it)->getNum());
+ };
+ #endif //_DEBUG
+ }
- nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- int dim = _solver->getMicroBC()->getDim();
- if (pbc != NULL && bgroup.size()>0){
- if (dim ==2){
- if (bgroup.size() == 4){
- groupIntersection(bgroup[0],bgroup[1],_bVertices);
- groupIntersection(bgroup[1],bgroup[2],_bVertices);
- groupIntersection(bgroup[2],bgroup[3],_bVertices);
- groupIntersection(bgroup[3],bgroup[0],_bVertices);
-
- l12.insert(bgroup[0]->vbegin(),bgroup[0]->vend());
- l41.insert(bgroup[1]->vbegin(),bgroup[1]->vend());
- l34.insert(bgroup[2]->vbegin(),bgroup[2]->vend());
- l23.insert(bgroup[3]->vbegin(),bgroup[3]->vend());
-
- #ifdef _DEBUG
- for (std::set<MVertex*>::iterator it = l12.begin(); it!= l12.end(); it++){
- Msg::Info("vertex %d 12",(*it)->getNum());
- };
- for (std::set<MVertex*>::iterator it = l23.begin(); it!= l23.end(); it++){
- Msg::Info("vertex %d 23",(*it)->getNum());
- };
- for (std::set<MVertex*>::iterator it = l34.begin(); it!= l34.end(); it++){
- Msg::Info("vertex %d 34",(*it)->getNum());
- };
- for (std::set<MVertex*>::iterator it = l41.begin(); it!= l41.end(); it++){
- Msg::Info("vertex %d 41",(*it)->getNum());
- };
- #endif //_DEBUG
-
- for (std::set<MVertex*>::iterator it = _bVertices.begin(); it!= _bVertices.end(); it++){
- _cVertices.insert(*it);
- }
- }
- else if (bgroup.size() == 2){
- l12.insert(bgroup[0]->vbegin(),bgroup[0]->vend());
- l34.insert(bgroup[2]->vbegin(),bgroup[2]->vend());
-
- #ifdef _DEBUG
- for (std::set<MVertex*>::iterator it = l12.begin(); it!= l12.end(); it++){
- Msg::Info("vertex %d 12",(*it)->getNum());
- };
- for (std::set<MVertex*>::iterator it = l34.begin(); it!= l34.end(); it++){
- Msg::Info("vertex %d 34",(*it)->getNum());
- };
- #endif //_DEBUG
- }
- else{
- Msg::Fatal("Error in setting PBC, the number of physical for PBC must be equal to 2, 4. In this case: %d ",bgroup.size());
- }
- }
- else if (dim ==3){
- if (bgroup.size() == 6){
- f1234 = bgroup[0];
- f5678 = bgroup[3];
- f1562 = bgroup[1];
- f4873 = bgroup[4];
- f1584 = bgroup[2];
- f2673 = bgroup[5];
-
- groupIntersection(f1234,f1562,l12);
- for (std::set<MVertex*>::iterator it = l12.begin(); it!= l12.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 12",(*it)->getNum());
- #endif //_DEBUG
- };
- groupIntersection(f1234,f2673,l23);
- for (std::set<MVertex*>::iterator it = l23.begin(); it!= l23.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 23",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f1234,f4873,l34);
- for (std::set<MVertex*>::iterator it = l34.begin(); it!= l34.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 34",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f1234,f1584,l41);
- for (std::set<MVertex*>::iterator it = l41.begin(); it!= l41.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 41",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f5678,f1562,l56);
- for (std::set<MVertex*>::iterator it = l56.begin(); it!= l56.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 56",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f5678,f2673,l67);
- for (std::set<MVertex*>::iterator it = l67.begin(); it!= l67.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 67",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f5678,f4873,l78);
- for (std::set<MVertex*>::iterator it = l78.begin(); it!= l78.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 78",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f5678,f1584,l85);
- for (std::set<MVertex*>::iterator it = l85.begin(); it!= l85.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 85",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f1562,f1584,l15);
- for (std::set<MVertex*>::iterator it = l15.begin(); it!= l15.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 15",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f1562,f2673,l26);
- for (std::set<MVertex*>::iterator it = l26.begin(); it!= l26.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 26",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f2673,f4873,l37);
- for (std::set<MVertex*>::iterator it = l37.begin(); it!= l37.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 37",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f4873,f1584,l48);
- for (std::set<MVertex*>::iterator it = l48.begin(); it!= l48.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 48",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(l12,l23,_cVertices);
- groupIntersection(l23,l34,_cVertices);
- groupIntersection(l34,l41,_cVertices);
- groupIntersection(l41,l12,_cVertices);
- groupIntersection(l56,l67,_cVertices);
- groupIntersection(l67,l78,_cVertices);
- groupIntersection(l78,l85,_cVertices);
- groupIntersection(l85,l56,_cVertices);
- }
- else if (bgroup.size() == 4){
- f1234 = bgroup[0];
- f5678 = bgroup[2];
-
- f1584 = bgroup[1];
- f2673 = bgroup[3];
-
- groupIntersection(f1234,f1584,l41);
- for (std::set<MVertex*>::iterator it = l41.begin(); it!= l41.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 41",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f1234,f2673,l23);
- for (std::set<MVertex*>::iterator it = l23.begin(); it!= l23.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 23",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f5678,f2673,l67);
- for (std::set<MVertex*>::iterator it = l67.begin(); it!= l67.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 67",(*it)->getNum());
- #endif //_DEBUG
- };
-
- groupIntersection(f5678,f1584,l85);
- for (std::set<MVertex*>::iterator it = l85.begin(); it!= l85.end(); it++){
- _bVertices.insert(*it);
- #ifdef _DEBUG
- Msg::Info("vertex %d 85",(*it)->getNum());
- #endif //_DEBUG
- };
-
- }
- else if (bgroup.size() == 2){
- f1234 = bgroup[0];
- f5678 = bgroup[1];
- }
- else{
- Msg::Fatal("Error in setting PBC, the number of physical for PBC must be equal to 2, 4, 6. In this case: %d ",bgroup.size());
- }
- }
-
- #ifdef _DEBUG
- for (std::set<MVertex*>::iterator it = _bVertices.begin(); it!= _bVertices.end(); it++){
- MVertex* v = *it;
- Msg::Info("Boundary vertex = %d",v->getNum());
- }
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- MVertex* v = *it;
- Msg::Info("Corner vertex = %d",v->getNum());
- }
- #endif //_DEBUG
-
- addControlPoints();
- addOtherCornerPoints();
- }
-};
-
-
-void pbcConstraintElementGroup::addControlPoints(){
+ // add control points, v1-v2-v4 in 2D, v1, v2, v4, v5 in 3D
std::set<MVertex*> vv;
- groupIntersection(l12,l41,vv);
- if (_solver->getMicroBC()->getDim() == 3){
- groupIntersection(l12,l15,vv);
- groupIntersection(l41,l15,vv);
- }
- if (vv.size()==0){
- MVertex* v1(NULL),*v2(NULL),*v3(NULL);
- if (l12.size()>0 and l41.size()>2 ){
- v1 = *(l12.begin());
- std::set<MVertex*>::iterator it = l41.begin();
- v2 = *(it);
- it++;
- v3 = *(it);
- }
- if (v1!=NULL and v2!=NULL and v3!=NULL){
- SPoint3 pp= project(v1->point(),v2->point(),v3->point());
+ if (dim ==2){
+ vv.clear();
+ InterpolationOperations::groupIntersection(l12,l41,vv);
+ if (vv.size() > 0){
+ _v1 = *(vv.begin());
+ }
+ else{
+ MElement* ex = *(bgroup[0]->begin());
+ MElement* ey = *(bgroup[1]->begin());
+ std::vector<MVertex*> vx, vy;
+ ex->getVertices(vx);
+ ey->getVertices(vy);
+ SPoint3 pp = computeTwoLineIntersectionOnSurface(vx[0]->point(),vx[1]->point(),vy[0]->point(),vy[1]->point());
_v1 = new MVertex(pp.x(),pp.y(),pp.z());
Msg::Info("Create root vertex v1 %d: %e %e %e",_v1->getNum(),pp.x(),pp.y(),pp.z());
- l12.insert(_v1);
- l41.insert(_v1);
- if (_solver->getMicroBC()->getDim()==3){
- l15.insert(_v1);
- }
_virtualVertices.insert(_v1);
_cVertices.insert(_v1);
_bVertices.insert(_v1);
}
- }
- else {
- _v1 = *(vv.begin());
- }
-
- vv.clear();
- groupIntersection(l12,l23,vv);
- if (_solver->getMicroBC()->getDim()==3){
- groupIntersection(l12,l26,vv);
- groupIntersection(l23,l26,vv);
- }
- if (vv.size()==0){
- MVertex* v1(NULL),*v2(NULL),*v3(NULL);
- if (l12.size()>0 and l23.size()>2){
- v1 = *(l12.begin());
- std::set<MVertex*>::iterator it = l23.begin();
- v2 = *(it);
- it++;
- v3 = *(it);
- }
- if (v1!=NULL and v2!=NULL and v3!=NULL){
- SPoint3 pp= project(v1->point(),v2->point(),v3->point());
+
+ // find v2
+ vv.clear();
+ InterpolationOperations::groupIntersection(l12,l23,vv);
+ if (vv.size() > 0){
+ _v2 = *(vv.begin());
+ }
+ else{
+ MElement* ex = *(bgroup[0]->begin());
+ MElement* ey = *(bgroup[3]->begin());
+ std::vector<MVertex*> vx, vy;
+ ex->getVertices(vx);
+ ey->getVertices(vy);
+ SPoint3 pp = computeTwoLineIntersectionOnSurface(vx[0]->point(),vx[1]->point(),vy[0]->point(),vy[1]->point());
_v2 = new MVertex(pp.x(),pp.y(),pp.z());
Msg::Info("Create root vertex v2 %d: %e %e %e",_v2->getNum(),pp.x(),pp.y(),pp.z());
- l12.insert(_v2);
_virtualVertices.insert(_v2);
_cVertices.insert(_v2);
_bVertices.insert(_v2);
}
- }
- else {
- _v2 = *(vv.begin());
- }
-
- vv.clear();
- groupIntersection(l34,l41,vv);
- if (_solver->getMicroBC()->getDim()==3){
- groupIntersection(l34,l48,vv);
- groupIntersection(l41,l48,vv);
- }
- if (vv.size()==0){
- MVertex* v1(NULL),*v2(NULL),*v3(NULL);
- if (l41.size()>0 and l34.size()>2){
- v1 = *(l41.begin());
- std::set<MVertex*>::iterator it = l34.begin();
- v2 = *(it);
- it++;
- v3 = *(it);
- }
- if (v1!=NULL and v2!=NULL and v3!=NULL){
- SPoint3 pp= project(v1->point(),v2->point(),v3->point());
+
+ // find v3
+ vv.clear();
+ InterpolationOperations::groupIntersection(l23,l34,vv);
+ if (vv.size()>0){
+ _v3 = *(vv.begin());
+ }
+ else{
+ MElement* ex = *(bgroup[3]->begin());
+ MElement* ey = *(bgroup[2]->begin());
+ std::vector<MVertex*> vx, vy;
+ ex->getVertices(vx);
+ ey->getVertices(vy);
+ SPoint3 pp = computeTwoLineIntersectionOnSurface(vx[0]->point(),vx[1]->point(),vy[0]->point(),vy[1]->point());
+ _v3 = new MVertex(pp.x(),pp.y(),pp.z());
+ Msg::Info("Create root vertex v3 %d: %e %e %e",_v3->getNum(),pp.x(),pp.y(),pp.z());
+ _virtualVertices.insert(_v3);
+ _cVertices.insert(_v3);
+ _bVertices.insert(_v3);
+ }
+
+ // find v4
+ vv.clear();
+ InterpolationOperations::groupIntersection(l34,l41,vv);
+ if (vv.size() > 0){
+ _v4 = *(vv.begin());
+ }
+ else{
+ MElement* ex = *(bgroup[2]->begin());
+ MElement* ey = *(bgroup[1]->begin());
+ std::vector<MVertex*> vx, vy;
+ ex->getVertices(vx);
+ ey->getVertices(vy);
+ SPoint3 pp = computeTwoLineIntersectionOnSurface(vx[0]->point(),vx[1]->point(),vy[0]->point(),vy[1]->point());
_v4 = new MVertex(pp.x(),pp.y(),pp.z());
Msg::Info("Create root vertex v4 %d: %e %e %e",_v4->getNum(),pp.x(),pp.y(),pp.z());
- l41.insert(_v4);
_virtualVertices.insert(_v4);
_cVertices.insert(_v4);
_bVertices.insert(_v4);
}
+
+ _periodic12[0] = _v2->x() - _v1->x();
+ _periodic12[1] = _v2->y() - _v1->y();
+ _periodic12[2] = _v2->z() - _v1->z();
+
+ _periodic14[0] = _v4->x() - _v1->x();
+ _periodic14[1] = _v4->y() - _v1->y();
+ _periodic14[2] = _v4->z() - _v1->z();
+
+ _periodic15 = crossprod(_periodic12,_periodic14);
+ _periodic15.normalize();
}
- else {
- _v4 = *(vv.begin());
- }
- if (_solver->getMicroBC()->getDim()==3){
+ else if (dim == 3){
vv.clear();
- groupIntersection(l85,l15,vv);
- groupIntersection(l85,l56,vv);
- groupIntersection(l56,l15,vv);
-
- if (vv.size()==0){
- MVertex* v1(NULL),*v2(NULL),*v3(NULL);
- if (l15.size()>0 and l85.size()>2){
- v1 = *(l15.begin());
- std::set<MVertex*>::iterator it = l85.begin();
- v2 = *(it);
- it++;
- v3 = *(it);
- }
-
- if (v1!=NULL and v2!=NULL and v3!=NULL){
- SPoint3 pp= project(v1->point(),v2->point(),v3->point());
- _v5 = new MVertex(pp.x(),pp.y(),pp.z());
- Msg::Info("Create root vertex v5 %d: %e %e %e",_v5->getNum(),pp.x(),pp.y(),pp.z());
- l15.insert(_v5);
- _virtualVertices.insert(_v5);
- _cVertices.insert(_v5);
- _bVertices.insert(_v5);
- }
-
+ InterpolationOperations::groupIntersection(l12,l41,vv);
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l12,l15,vv);
}
- else {
- _v5 = *(vv.begin());
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l41,l15,vv);
}
-
- if (_v1==NULL){
+ if (vv.size() > 0){
+ _v1 = *(vv.begin());
+ }
+ else{
MElement* e1 = *(f1234->begin());
MElement* e2 = *(f1584->begin());
MElement* e3 = *(f1562->begin());
@@ -608,8 +801,20 @@ void pbcConstraintElementGroup::addControlPoints(){
_cVertices.insert(_v1);
_bVertices.insert(_v1);
}
-
- if (_v2==NULL){
+
+ // v2
+ vv.clear();
+ InterpolationOperations::groupIntersection(l12,l23,vv);
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l12,l26,vv);
+ }
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l23,l26,vv);
+ }
+ if (vv.size() > 0) {
+ _v2 = *(vv.begin());
+ }
+ else{
MElement* e1 = *(f1234->begin());
MElement* e2 = *(f2673->begin());
MElement* e3 = *(f1562->begin());
@@ -636,8 +841,20 @@ void pbcConstraintElementGroup::addControlPoints(){
_cVertices.insert(_v2);
_bVertices.insert(_v2);
}
-
- if (_v4==NULL){
+
+ // find v4
+ vv.clear();
+ InterpolationOperations::groupIntersection(l34,l41,vv);
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l34,l48,vv);
+ }
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l41,l48,vv);
+ }
+ if (vv.size() > 0) {
+ _v4 = *(vv.begin());
+ }
+ else{
MElement* e1 = *(f1234->begin());
MElement* e2 = *(f1584->begin());
MElement* e3 = *(f4873->begin());
@@ -664,8 +881,20 @@ void pbcConstraintElementGroup::addControlPoints(){
_cVertices.insert(_v4);
_bVertices.insert(_v4);
}
-
- if (_v5==NULL){
+
+ // find v5
+ vv.clear();
+ InterpolationOperations::groupIntersection(l85,l15,vv);
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l85,l56,vv);
+ }
+ if (vv.size() == 0){
+ InterpolationOperations::groupIntersection(l56,l15,vv);
+ }
+ if (vv.size() >0) {
+ _v5 = *(vv.begin());
+ }
+ else{
MElement* e1 = *(f1584->begin());
MElement* e2 = *(f5678->begin());
MElement* e3 = *(f1562->begin());
@@ -692,72 +921,39 @@ void pbcConstraintElementGroup::addControlPoints(){
_cVertices.insert(_v5);
_bVertices.insert(_v5);
}
-
- }
-
- _periodic12[0] = _v2->x() - _v1->x();
- _periodic12[1] = _v2->y() - _v1->y();
- _periodic12[2] = _v2->z() - _v1->z();
-
- _periodic14[0] = _v4->x() - _v1->x();
- _periodic14[1] = _v4->y() - _v1->y();
- _periodic14[2] = _v4->z() - _v1->z();
-
- if (_solver->getMicroBC()->getDim() == 3){
- _periodic15[0] = _v5->x() - _v1->x();
- _periodic15[1] = _v5->y() - _v1->y();
- _periodic15[2] = _v5->z() - _v1->z();
- }
-
- _periodic12.print("periodic 1 --> 2");
- _periodic14.print("periodic 1 --> 4");
-
-
- if (_solver->getMicroBC()->getDim() == 3){
- _periodic15.print("periodic 1 --> 5");
- }
-
-};
-
-void pbcConstraintElementGroup::addOtherCornerPoints(){
- std::set<MVertex*> vv;
- // find _v3
- vv.clear();
- groupIntersection(l23,l34,vv);
- if (_solver->getMicroBC()->getDim() ==3){
- if (vv.size()==0)
- groupIntersection(l34,l37,vv);
- if (vv.size() ==0)
- groupIntersection(l23,l37,vv);
- }
-
- if (vv.size()==0){
- double x = _v4->x() + _v2->x() - _v1->x();
- double y = _v4->y() + _v2->y() - _v1->y();
- double z = _v4->z() + _v2->z() - _v1->z();
- _v3 = new MVertex(x,y,z);
- Msg::Info("Create root vertex v3 %d: %e %e %e",_v3->getNum(),x,y,z);
- l23.insert(_v3);
- l34.insert(_v3);
- if (_solver->getMicroBC()->getDim() ==3)
- l37.insert(_v3);
- _virtualVertices.insert(_v3);
- _cVertices.insert(_v3);
- _bVertices.insert(_v3);
- }
- else {
- _v3 = *(vv.begin());
- }
-
- if (_solver->getMicroBC()->getDim() ==3){
+
+
+ // find v3
+ vv.clear();
+ InterpolationOperations::groupIntersection(l23,l34,vv);
+ if (vv.size()==0){
+ InterpolationOperations::groupIntersection(l34,l37,vv);
+ }
+ if (vv.size() ==0){
+ InterpolationOperations::groupIntersection(l23,l37,vv);
+ }
+ if (vv.size() > 0){
+ _v3 = *(vv.begin());
+ }
+ else{
+ double x = _v4->x() + _v2->x() - _v1->x();
+ double y = _v4->y() + _v2->y() - _v1->y();
+ double z = _v4->z() + _v2->z() - _v1->z();
+ _v3 = new MVertex(x,y,z);
+ Msg::Info("Create root vertex v3 %d: %e %e %e",_v3->getNum(),x,y,z);
+ _virtualVertices.insert(_v3);
+ _cVertices.insert(_v3);
+ _bVertices.insert(_v3);
+ }
+
// find _v6
vv.clear();
- groupIntersection(l26,l56,vv);
+ InterpolationOperations::groupIntersection(l26,l56,vv);
if (vv.size() == 0)
- groupIntersection(l56,l67,vv);
+ InterpolationOperations::groupIntersection(l56,l67,vv);
if (vv.size() == 0)
- groupIntersection(l26,l67,vv);
+ InterpolationOperations::groupIntersection(l26,l67,vv);
if (vv.size()==0){
double x = _v2->x() + _v5->x() - _v1->x();
@@ -765,9 +961,6 @@ void pbcConstraintElementGroup::addOtherCornerPoints(){
double z = _v2->z() + _v5->z() - _v1->z();
_v6 = new MVertex(x,y,z);
Msg::Info("Create root vertex v6 %d: %e %e %e",_v6->getNum(),x,y,z);
- l26.insert(_v6);
- l56.insert(_v6);
- l67.insert(_v6);
_virtualVertices.insert(_v6);
_cVertices.insert(_v6);
_bVertices.insert(_v6);
@@ -778,12 +971,12 @@ void pbcConstraintElementGroup::addOtherCornerPoints(){
// find _v8
vv.clear();
- groupIntersection(l85,l48,vv);
+ InterpolationOperations::groupIntersection(l85,l48,vv);
if (vv.size() == 0)
- groupIntersection(l48,l78,vv);
+ InterpolationOperations::groupIntersection(l48,l78,vv);
if (vv.size() == 0)
- groupIntersection(l85,l78,vv);
+ InterpolationOperations::groupIntersection(l85,l78,vv);
if (vv.size()==0){
double x = _v4->x() + _v5->x() - _v1->x();
@@ -791,9 +984,6 @@ void pbcConstraintElementGroup::addOtherCornerPoints(){
double z = _v4->z() + _v5->z() - _v1->z();
_v8 = new MVertex(x,y,z);
Msg::Info("Create root vertex v8 %d: %e %e %e",_v8->getNum(),x,y,z);
- l48.insert(_v8);
- l78.insert(_v8);
- l85.insert(_v8);
_virtualVertices.insert(_v8);
_cVertices.insert(_v8);
_bVertices.insert(_v8);
@@ -804,12 +994,12 @@ void pbcConstraintElementGroup::addOtherCornerPoints(){
// find _v7
vv.clear();
- groupIntersection(l78,l37,vv);
+ InterpolationOperations::groupIntersection(l78,l37,vv);
if (vv.size() == 0)
- groupIntersection(l37,l67,vv);
+ InterpolationOperations::groupIntersection(l37,l67,vv);
if (vv.size() == 0)
- groupIntersection(l78,l67,vv);
+ InterpolationOperations::groupIntersection(l78,l67,vv);
if (vv.size()==0){
double x = _v3->x() + _v5->x() - _v1->x();
@@ -817,9 +1007,6 @@ void pbcConstraintElementGroup::addOtherCornerPoints(){
double z = _v3->z() + _v5->z() - _v1->z();
_v7 = new MVertex(x,y,z);
Msg::Info("Create root vertex v7 %d: %e %e %e",_v7->getNum(),x,y,z);
- l78.insert(_v7);
- l37.insert(_v7);
- l67.insert(_v7);
_virtualVertices.insert(_v7);
_cVertices.insert(_v7);
_bVertices.insert(_v7);
@@ -827,333 +1014,26 @@ void pbcConstraintElementGroup::addOtherCornerPoints(){
else {
_v7 = *(vv.begin());
}
- }
-};
-
-void pbcConstraintElementGroup::createPolynomialBasePoints(){
- int dim = _solver->getMicroBC()->getDim();
- nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- int degreeX = pbc->getPBCPolynomialDegree(0);
- int degreeY = pbc->getPBCPolynomialDegree(1);;
- int degreeZ = pbc->getPBCPolynomialDegree(2);;
-
- if (!pbc->isForcedInterpolationDegree()){
- int sizeX = l12.size();
- int sizeY = l41.size();
- int sizeZ = l15.size();
-
- Msg::Info("Maximal number of Boundary Nodes in X direction %d, Y direction %d, Z direction %d",sizeX,sizeY,sizeZ);
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::LIM or
- pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN){
- // first order
- if (degreeX>sizeX-1) degreeX = sizeX-1;
- if (degreeY>sizeY-1) degreeY = sizeY-1;
- if (degreeZ>sizeZ-1 && dim == 3) degreeZ = sizeZ-1;
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM or
- pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
- // second-order
- if (degreeX>(sizeX-1)/2) degreeX = (sizeX-1)/2;
- if (degreeY>(sizeY-1)/2) degreeY = (sizeY-1)/2;
- if (degreeZ>(sizeZ-1)/2 && dim == 3) degreeZ = (sizeZ-1)/2;
- }
- else{
- Msg::Fatal("interpolation type must be defined");
- }
- if (degreeX<1) degreeX = 1;
- if (degreeY<1) degreeY = 1;
- if (degreeZ<1) degreeZ = 1;
-
- // modify degree
- pbc->setPBCPolynomialDegree(0,degreeX);
- pbc->setPBCPolynomialDegree(1,degreeY);
- pbc->setPBCPolynomialDegree(2,degreeZ);
- }
-
- Msg::Info("Interpolation degree used in X direction %d, Y direction %d, Z direction %d",degreeX,degreeY,degreeZ);
-
- bool ok = false;
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN or
- pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
- ok = true;
- }
- if (pbc->addNewVertices(0) or l12.size()<=2 or ok){
- this->createNewVerticesForInterpolationPoints(degreeX,_v1,_v2,_xBase,_cVertices);
- }
- else
- this->getInterpolatioPointsFromExistingVertices(degreeX,_v1,_v4,l12.begin(), l12.end(),_xBase);
-
-
- if (pbc->addNewVertices(1) or l41.size()<=2 or ok){
- this->createNewVerticesForInterpolationPoints(degreeY,_v1,_v4,_yBase,_cVertices);
- }
- else
- this->getInterpolatioPointsFromExistingVertices(degreeY,_v1,_v2,l41.begin(), l41.end(),_yBase);
-
- if (dim == 3){
- if (pbc->addNewVertices(2) or l15.size()<=2 or ok){
- this->createNewVerticesForInterpolationPoints(degreeZ,_v1,_v5,_zBase,_cVertices);
- }
- else
- this->getInterpolatioPointsFromExistingVertices(degreeZ,_v1,_v4,l15.begin(), l15.end(),_zBase);
- }
-};
-
-void pbcConstraintElementGroup::createBCGroup3D( std::vector<MVertex*>& vlx, std::vector<MVertex*>& vly,
- MVertex* vc, groupOfElements* g){
- if (vlx[0] != vly[0]){
- Msg::Fatal("two BC group must be started by same vertex");
- return;
- }
- MVertex* vroot = vlx[0];
- nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- if (pbc == NULL) return;
- int NX = vlx.size();
- int NY = vly.size();
-
- std::vector<std::vector<MVertex*> > vmat;
- vmat.push_back(vlx);
-
- for (int j=1; j<NY; j++){
- std::vector<MVertex*> vv;
- vv.push_back(vly[j]);
- for (int i=1; i< NX; i++){
- MVertex* vertex = NULL;
- if (j == NY-1 && i == NX-1){
- vertex = vc;
- }
- else{
- vertex = new MVertex(0,0,0);
- vertex->x() = vlx[i]->x() - vroot->x()+ vly[j]->x();
- vertex->y() = vlx[i]->y() - vroot->y()+ vly[j]->y();
- vertex->z() = vlx[i]->z() - vroot->z()+ vly[j]->z();
- _virtualVertices.insert(vertex);
- }
- vv.push_back(vertex);
- }
- vmat.push_back(vv);
- }
-
-
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN){
- for (int i=0; i<NX-1; i++){
- for (int j=0; j<NY-1; j++){
- MElement* e = new MQuadrangle(vmat[j][i],vmat[j][i+1],vmat[j+1][i+1],vmat[j+1][i]);
- g->insert(e);
- }
- }
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
- for (int i=0; i<NX-1; i+=2){
- for (int j=0; j<NY-1; j+=2){
- MElement* e = new MQuadrangle9(vmat[j][i],vmat[j][i+2],vmat[j+2][i+2],vmat[j+2][i],
- vmat[j][i+1],vmat[j+1][i+2],vmat[j+2][i+1],vmat[j+1][i],
- vmat[j+1][i+1]);
- g->insert(e);
- }
- }
- }
-};
-
-void pbcConstraintElementGroup::intersectionBetweenTwoFaces(const groupOfElements* face1, const groupOfElements* face2, groupOfElements* & line) const{
-
- std::vector<std::vector<MVertex*> > allEdges;
- for (groupOfElements::elementContainer::iterator itL = face1->begin(); itL!= face1->end(); itL++){
- MElement* eleL = *itL;
- if (eleL->getDim() != 2){
- Msg::Fatal("This fuction is implemented from 2D elements only");
- }
- bool found = false;
- for (int i=0; i< eleL->getNumEdges(); i++){
- std::vector<MVertex*> vvL;
- eleL->getEdgeVertices(i,vvL);
- // found in face3
- for (groupOfElements::elementContainer::iterator itR = face2->begin(); itR!= face2->end(); itR++){
- MElement* eleR = *itR;
- for (int j=0; j< eleR->getNumEdges(); j++){
- std::vector<MVertex*> vvR;
- eleR->getEdgeVertices(j,vvR);
- if (((vvL[0]->getNum() == vvR[0]->getNum()) and (vvL[1]->getNum() == vvR[1]->getNum())) or
- ((vvL[0]->getNum() == vvR[1]->getNum()) and (vvL[1]->getNum() == vvR[0]->getNum()))){
- allEdges.push_back(vvL);
- found = true;
- break;
- }
-
- }
- if (found) break;
- }
- if (found) break;
-
- }
- }
-
- // add to grEle
- if (line == NULL) line = new groupOfElements();
- for (int i=0; i< allEdges.size(); i++){
- std::vector<MVertex*>& vv = allEdges[i];
- MElement* ele = new MLineN(vv);
- line->insert(ele);
- }
-};
-
-void pbcConstraintElementGroup::intersectionBetweenTwoLines(const groupOfElements* line1, const groupOfElements* line2, groupOfElements* & pt) const{
- std::vector<MVertex*> allVertex;
- for (groupOfElements::vertexContainer::iterator itL = line1->vbegin(); itL!= line1->vend(); itL++){
- MVertex* v1 = *itL;
- for (groupOfElements::vertexContainer::iterator itR = line2->vbegin(); itR!= line2->vend(); itR++){
- MVertex* v2 = *itR;
- if (v1->getNum() == v2->getNum()){
- allVertex.push_back(v1);
- break;
- }
- }
- }
-
- if (pt == NULL) pt = new groupOfElements();
- for (int i=0; i< allVertex.size(); i++){
- MElement* ele = new MPoint(allVertex[i]);
- pt->insert(ele);
- }
-};
-
-void pbcConstraintElementGroup::createAllBCGroups(){
- int dim = _solver->getMicroBC()->getDim();
- if (_gX != NULL) delete _gX;
- _gX = new groupOfElements();
- if (_gY != NULL) delete _gY;
- _gY = new groupOfElements();
-
- nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN){
- for (int i=0; i<_xBase.size()-1; i++){
- MElement* line = new MLine(_xBase[i],_xBase[i+1]);
- _gX->insert(line);
- }
-
- for (int i=0; i<_yBase.size()-1; i++){
- MElement* line = new MLine(_yBase[i],_yBase[i+1]);
- _gY->insert(line);
- }
-
- if (dim == 3){
- _gZ = new groupOfElements();
- for (int i=0; i<_zBase.size()-1; i++){
- MElement* line = new MLine(_zBase[i],_zBase[i+1]);
- _gZ->insert(line);
- }
- }
- }
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
- std::vector<MVertex*> xtemp(_xBase);
- _xBase.clear();
- for (int i=0; i<xtemp.size()-1; i++){
- _xBase.push_back(xtemp[i]);
- MVertex* v = new MVertex(0.5*(xtemp[i+1]->x()+xtemp[i]->x()),
- 0.5*(xtemp[i+1]->y()+xtemp[i]->y()),
- 0.5*(xtemp[i+1]->z()+xtemp[i]->z()));
-
- _virtualVertices.insert(v);
- _bVertices.insert(v);
- _xBase.push_back(v);
-
- MElement* line = new MLine3(xtemp[i],xtemp[i+1],v);
- _gX->insert(line);
- if (i == xtemp.size()-2){
- _xBase.push_back(xtemp[i+1]);
- }
- }
-
- std::vector<MVertex*> ytemp(_yBase);
- _yBase.clear();
- for (int i=0; i<ytemp.size()-1; i++){
- _yBase.push_back(ytemp[i]);
- MVertex* v = new MVertex(0.5*(ytemp[i+1]->x()+ytemp[i]->x()),
- 0.5*(ytemp[i+1]->y()+ytemp[i]->y()),
- 0.5*(ytemp[i+1]->z()+ytemp[i]->z()));
-
- _virtualVertices.insert(v);
- _bVertices.insert(v);
- _yBase.push_back(v);
-
- MElement* line = new MLine3(ytemp[i],ytemp[i+1],v);
- _gY->insert(line);
- if (i == ytemp.size()-2){
- _yBase.push_back(ytemp[i+1]);
- }
- }
-
- if (dim== 3){
- std::vector<MVertex*> ztemp(_zBase);
- _zBase.clear();
- if (_gZ != NULL) delete _gZ;
- _gZ = new groupOfElements();
- for (int i=0; i<ztemp.size()-1; i++){
- _zBase.push_back(ztemp[i]);
- MVertex* v = new MVertex(0.5*(ztemp[i+1]->x()+ztemp[i]->x()),
- 0.5*(ztemp[i+1]->y()+ztemp[i]->y()),
- 0.5*(ztemp[i+1]->z()+ztemp[i]->z()));
-
- _virtualVertices.insert(v);
- _bVertices.insert(v);
- _zBase.push_back(v);
-
- MElement* line = new MLine3(ztemp[i],ztemp[i+1],v);
- _gZ->insert(line);
- if (i == ztemp.size()-2){
- _zBase.push_back(ztemp[i+1]);
- }
- }
+
+
+ _periodic12[0] = _v2->x() - _v1->x();
+ _periodic12[1] = _v2->y() - _v1->y();
+ _periodic12[2] = _v2->z() - _v1->z();
- }
- };
+ _periodic14[0] = _v4->x() - _v1->x();
+ _periodic14[1] = _v4->y() - _v1->y();
+ _periodic14[2] = _v4->z() - _v1->z();
- if (dim == 3){
- if (_gXY != NULL) delete _gXY;
- if (_gYZ != NULL) delete _gYZ;
- if (_gZX != NULL) delete _gZX;
- _gXY = new groupOfElements();
- _gZX = new groupOfElements();
- _gYZ = new groupOfElements();
-
- this->createBCGroup3D(_xBase,_yBase,_v3,_gXY);
- this->createBCGroup3D(_yBase,_zBase,_v8,_gYZ);
- this->createBCGroup3D(_zBase,_xBase,_v6,_gZX);
+ _periodic15[0] = _v5->x() - _v1->x();
+ _periodic15[1] = _v5->y() - _v1->y();
+ _periodic15[2] = _v5->z() - _v1->z();
}
+ _periodic12.print("periodic 1 --> 2");
+ _periodic14.print("periodic 1 --> 4");
+ _periodic15.print("periodic 1 --> 5");
+
};
-template<class Iterator>
-MVertex* pbcConstraintElementGroup::findNearestVertex(Iterator itbegin, Iterator itend, MVertex* vp){
- MVertex* vinit = *itbegin;
- double dist=vinit->distance(vp);
- for (Iterator it=itbegin; it!=itend; it++){
- MVertex* v=*it;
- double tempDist=v->distance(vp);
- if (tempDist<dist) {
- dist=tempDist;
- vinit=v;
- };
- };
- return vinit;
-};
-
-MVertex* pbcConstraintElementGroup::findNearestVertex(groupOfElements* g, MVertex* vp){
- /*
- groupOfElements::vertexContainer::iterator it = g->vbegin();
- MVertex* vinit = *it;
- double dist=vinit->distance(vp);
- for (groupOfElements::vertexContainer::iterator it=g->vbegin(); it!=g->vend(); it++){
- MVertex* v=*it;
- double tempDist=v->distance(vp);
- if (tempDist<dist) {
- dist=tempDist;
- vinit=v;
- };
- };
- return vinit;
- */
- return findNearestVertex(g->vbegin(), g->vend(),vp);
-};
template<class Iterator>
void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(Iterator posivtiveitbegin, Iterator positiveitend,
@@ -1162,7 +1042,7 @@ void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(Iterator posivt
for (Iterator itp = posivtiveitbegin; itp!= positiveitend; itp++){
MVertex* v = *itp;
if ( others.find(v) == others.end()){
- MVertex* vn = findNearestVertex(negativeitbegin,negativeitend,v);
+ MVertex* vn = InterpolationOperations::findNearestVertex(negativeitbegin,negativeitend,v);
constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,-1,v,vn,vrootP,vrootN);
_allConstraint.insert(cel);
};
@@ -1177,7 +1057,7 @@ void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(const int comp,
for (Iterator itp = posivtiveitbegin; itp!= positiveitend; itp++){
MVertex* v = *itp;
if ( others.find(v) == others.end()){
- MVertex* vn = findNearestVertex(negativeitbegin,negativeitend,v);
+ MVertex* vn = InterpolationOperations::findNearestVertex(negativeitbegin,negativeitend,v);
constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,comp,v,vn,vrootP,vrootN,fact);
_allConstraint.insert(cel);
};
@@ -1185,77 +1065,6 @@ void pbcConstraintElementGroup::periodicConditionForPeriodicMesh(const int comp,
};
-
-void pbcConstraintElementGroup::createNewVerticesForInterpolationPoints(int degree, MVertex* vinit, MVertex* vend,
- std::vector<MVertex*>& base,
- std::set<MVertex*>& others){
- base.clear();
- double x0 = vinit->x(); double x1 = vend->x();
- double y0 = vinit->y(); double y1 = vend->y();
- double z0 = vinit->z(); double z1 = vend->z();
- if (others.find(vinit) == others.end()){
- MVertex* v = new MVertex(x0,y0,z0);
- base.push_back(v);
- _virtualVertices.insert(v);
- }
- else{
- base.push_back(vinit);
- };
- double dx = (x1-x0)/degree;
- double dy = (y1-y0)/degree;
- double dz = (z1-z0)/degree;
- for (int j=0; j<degree-1; j++){
- double xi = x0+ (j+1)*dx;
- double yi = y0+ (j+1)*dy;
- double zi = z0+ (j+1)*dz;
- MVertex* v = new MVertex(xi,yi,zi);
- base.push_back(v);
- _virtualVertices.insert(v);
- };
-
- if (others.find(vend) == others.end()){
- MVertex* v = new MVertex(x1,y1,z1);
- base.push_back(v);
- _virtualVertices.insert(v);
- }
- else{
- base.push_back(vend);
- };
-};
-
-template<class Iterator>
-void pbcConstraintElementGroup::getInterpolatioPointsFromExistingVertices(int degree, MVertex* v1, MVertex* v2,
- Iterator itbegin, Iterator itend,
- std::vector<MVertex*>& base){
- base.clear();
- std::vector<MVertex*> list;
- std::vector<double> vals;
- for (Iterator it=itbegin; it!=itend; it++){
- MVertex* v=*it;
- list.push_back(v);
- vals.push_back(distance(v->point(),v1->point(),v2->point()));
- };
-
- int size=list.size();
- for (int i=0; i<size; i++){
- for (int j=i; j<size; j++){
- if (vals[j]<vals[i]){
- double temp=vals[i]; vals[i]=vals[j]; vals[j]=temp;
- MVertex* vtemp=list[i]; list[i]=list[j]; list[j]=vtemp;
- };
- };
- };
- if (degree>size-1){
- degree=size-1;
- std::cout<<"Modify the imposed degree: "<<degree<<std::endl;
- };
- base.push_back(list[0]);
- for (int j=0; j<degree-1; j++){
- base.push_back(list[(j+1)*(size-1)/degree]);
- };
- base.push_back(list[size-1]);
-};
-
// for cubic spline
template<class Iterator>
void pbcConstraintElementGroup::cubicSplineFormCondition(Iterator itbegin, Iterator itend,
@@ -1586,30 +1395,35 @@ void pbcConstraintElementGroup::lagrangePeriodicConditionCoonsPatch(Iterator itb
};
};
-bool pbcConstraintElementGroup::isInside(MVertex* v, MElement* e) const{
- int dim = e->getDim();
- std::vector<MVertex*> vv;
- e->getVertices(vv);
-
- SPoint3 pp;
- if (dim == 1){
- pp = project(v->point(),vv[0]->point(), vv[1]->point());
- }
- else if (dim ==2){
- pp = project(v->point(),vv[0]->point(), vv[1]->point(),vv[2]->point());
- }
-
- double xyz[3] = {pp[0],pp[1],pp[2]};
- double uvw[3];
-
- e->xyz2uvw(xyz,uvw);
- if (e->isInside(uvw[0],uvw[1],uvw[2])){
- return true;
- }
- else return false;
-
+template<class Iterator>
+void pbcConstraintElementGroup::pbcFormFEConstraintElementGroup(Iterator itbegin, Iterator itend,
+ groupOfElements* g, std::set<MVertex*>& others,
+ MVertex* vrootP, MVertex* vrootN){
+ std::set<MVertex*> allCheck(others);
+ allCheck.insert(g->vbegin(),g->vend());
+ for (Iterator it = itbegin; it!= itend;it++){
+ MVertex* v = *it;
+ if (allCheck.find(v)==allCheck.end()){
+ MElement* e = NULL;
+ for (groupOfElements::elementContainer::iterator itg = g->begin(); itg!= g->end(); itg++){
+ MElement* ele = *itg;
+ if (InterpolationOperations::isInside(v,ele)){
+ e = ele;
+ break;
+ }
+ }
+ if (e == NULL){
+ Msg::Error("receive element not found");
+ }
+ else{
+ constraintElement* cele = new FEConstraintElement(_solver->getMicroBC() ,_LagSpace,_MultSpace, -1,v,e,vrootP,vrootN);
+ _allConstraint.insert(cele);
+ }
+ };
+ };
};
+
template<class Iterator>
void pbcConstraintElementGroup::pbcFEConstraintElementGroup(Iterator itbegin, Iterator itend,
groupOfElements* g, std::set<MVertex*>& others,
@@ -1620,7 +1434,7 @@ void pbcConstraintElementGroup::pbcFEConstraintElementGroup(Iterator itbegin, It
MElement* e = NULL;
for (groupOfElements::elementContainer::iterator itg = g->begin(); itg!= g->end(); itg++){
MElement* ele = *itg;
- if (isInside(v,ele)){
+ if (InterpolationOperations::isInside(v,ele)){
e = ele;
break;
}
@@ -1640,41 +1454,92 @@ void pbcConstraintElementGroup::createFEConstraintElementGroup(){
int dim = _solver->getMicroBC()->getDim();
nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- // create base
- this->createPolynomialBasePoints();
- this->createAllBCGroups();
-
+ // create base
+ std::set<MVertex*> newVirVertices;
+ InterpolationOperations::createPolynomialBasePoints(pbc,_v1,_v2,_v4,_v5,l12,l41,l15,_xBase,_yBase,_zBase,newVirVertices);
+
if (dim ==2){
- this->pbcFEConstraintElementGroup(l12.begin(),l12.end(),_gX,_cVertices,_v1,_v1);
+ if (_gX == NULL) _gX = new groupOfElements();
+ else _gX->clearAll();
+
+ if (_gY == NULL) _gY = new groupOfElements();
+ else _gY->clearAll();
+
+ if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN){
+ InterpolationOperations::createGroupLinearLineElementFromVerticesVector(_xBase,*_gX);
+ InterpolationOperations::createGroupLinearLineElementFromVerticesVector(_yBase,*_gY);
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
+ InterpolationOperations::createGroupQuadraticLineElementFromVerticesVector(_xBase,*_gX,newVirVertices);
+ InterpolationOperations::createGroupQuadraticLineElementFromVerticesVector(_yBase,*_gY,newVirVertices);
+ }
+ _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
+
+ this->pbcFormFEConstraintElementGroup(l12.begin(),l12.end(),_gX,_cVertices,_v1,_v1);
this->pbcFEConstraintElementGroup(l34.begin(),l34.end(),_gX,_cVertices,_v4,_v1);
- this->pbcFEConstraintElementGroup(l41.begin(),l41.end(),_gY,_cVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(l41.begin(),l41.end(),_gY,_cVertices,_v1,_v1);
this->pbcFEConstraintElementGroup(l23.begin(),l23.end(),_gY,_cVertices,_v2,_v1);
}
else if (dim == 3){
- #ifdef _DEBUG
- Msg::Error("ele num gxy =%d, gyz = %d, gzx = %d",_gXY->size(),_gYZ->size(), _gZX->size());
- #endif
-
+ if (_gX == NULL) _gX = new groupOfElements();
+ else _gX->clearAll();
+ if (_gY == NULL) _gY = new groupOfElements();
+ else _gY->clearAll();
+ if (_gZ == NULL) _gZ = new groupOfElements();
+ else _gZ->clearAll();
+ if (_gXY == NULL) _gXY = new groupOfElements();
+ else _gXY->clearAll();
+ if (_gYZ == NULL) _gYZ = new groupOfElements();
+ else _gYZ->clearAll();
+ if (_gZX == NULL) _gZX = new groupOfElements();
+ else _gZX->clearAll();
+
+
+ if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN){
+ InterpolationOperations::createGroupLinearLineElementFromVerticesVector(_xBase,*_gX);
+ InterpolationOperations::createGroupLinearLineElementFromVerticesVector(_yBase,*_gY);
+ InterpolationOperations::createGroupLinearLineElementFromVerticesVector(_zBase,*_gZ);
+
+ _bVertices.insert(newVirVertices.begin(),newVirVertices.end());
+
+ InterpolationOperations::createLinearElementBCGroup3DFromTwoEdges(_xBase,_yBase,_v3,*_gXY,newVirVertices);
+ InterpolationOperations::createLinearElementBCGroup3DFromTwoEdges(_yBase,_zBase,_v8,*_gYZ,newVirVertices);
+ InterpolationOperations::createLinearElementBCGroup3DFromTwoEdges(_zBase,_xBase,_v6,*_gZX,newVirVertices);
+ _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
+
+ }
+ else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
+ InterpolationOperations::createGroupQuadraticLineElementFromVerticesVector(_xBase,*_gX,newVirVertices);
+ InterpolationOperations::createGroupQuadraticLineElementFromVerticesVector(_yBase,*_gY,newVirVertices);
+ InterpolationOperations::createGroupQuadraticLineElementFromVerticesVector(_zBase,*_gZ,newVirVertices);
+ _bVertices.insert(newVirVertices.begin(),newVirVertices.end());
+
+ InterpolationOperations::createQuadraticElementBCGroup3DFromTwoEdges(_xBase,_yBase,_v3,*_gXY,newVirVertices);
+ InterpolationOperations::createQuadraticElementBCGroup3DFromTwoEdges(_yBase,_zBase,_v8,*_gYZ,newVirVertices);
+ InterpolationOperations::createQuadraticElementBCGroup3DFromTwoEdges(_zBase,_xBase,_v6,*_gZX,newVirVertices);
+ _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
+ }
+
- this->pbcFEConstraintElementGroup(l12.begin(),l12.end(),_gX,_cVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(l12.begin(),l12.end(),_gX,_cVertices,_v1,_v1);
this->pbcFEConstraintElementGroup(l34.begin(),l34.end(),_gX,_cVertices,_v4,_v1);
this->pbcFEConstraintElementGroup(l78.begin(),l78.end(),_gX,_cVertices,_v8,_v1);
this->pbcFEConstraintElementGroup(l56.begin(),l56.end(),_gX,_cVertices,_v5,_v1);
- this->pbcFEConstraintElementGroup(l41.begin(),l41.end(),_gY,_cVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(l41.begin(),l41.end(),_gY,_cVertices,_v1,_v1);
this->pbcFEConstraintElementGroup(l23.begin(),l23.end(),_gY,_cVertices,_v2,_v1);
this->pbcFEConstraintElementGroup(l67.begin(),l67.end(),_gY,_cVertices,_v6,_v1);
this->pbcFEConstraintElementGroup(l85.begin(),l85.end(),_gY,_cVertices,_v5,_v1);
- this->pbcFEConstraintElementGroup(l15.begin(),l15.end(),_gZ,_cVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(l15.begin(),l15.end(),_gZ,_cVertices,_v1,_v1);
this->pbcFEConstraintElementGroup(l26.begin(),l26.end(),_gZ,_cVertices,_v2,_v1);
this->pbcFEConstraintElementGroup(l37.begin(),l37.end(),_gZ,_cVertices,_v3,_v1);
this->pbcFEConstraintElementGroup(l48.begin(),l48.end(),_gZ,_cVertices,_v4,_v1);
- this->pbcFEConstraintElementGroup(f1234->vbegin(), f1234->vend(),_gXY,_bVertices,_v1,_v1);
- this->pbcFEConstraintElementGroup(f1584->vbegin(), f1584->vend(),_gYZ,_bVertices,_v1,_v1);
- this->pbcFEConstraintElementGroup(f1562->vbegin(), f1562->vend(),_gZX,_bVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(f1234->vbegin(), f1234->vend(),_gXY,_bVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(f1584->vbegin(), f1584->vend(),_gYZ,_bVertices,_v1,_v1);
+ this->pbcFormFEConstraintElementGroup(f1562->vbegin(), f1562->vend(),_gZX,_bVertices,_v1,_v1);
this->pbcFEConstraintElementGroup(f5678->vbegin(), f5678->vend(),_gXY,_bVertices,_v5,_v1);
this->pbcFEConstraintElementGroup(f2673->vbegin(), f2673->vend(),_gYZ,_bVertices,_v2,_v1);
@@ -1693,15 +1558,16 @@ void pbcConstraintElementGroup::createProjectConstraintElementGroup(){
this->pbcFEConstraintElementGroup(l23.begin(),l23.end(),bgroup[1],_cVertices,_v2,_v1);
}
else if (dim == 3){
- if (_gX) delete _gX;
- if (_gY) delete _gY;
- if (_gZ) delete _gZ;
- _gX = new groupOfElements();
- intersectionBetweenTwoFaces(f1234,f1562,_gX);
- _gY = new groupOfElements();
- intersectionBetweenTwoFaces(f1234,f1584,_gY);
- _gZ = new groupOfElements();
- intersectionBetweenTwoFaces(f1562,f1584,_gZ);
+ if (_gX == NULL) _gX = new groupOfElements();
+ else _gX->clearAll();
+ if (_gY == NULL) _gY = new groupOfElements();
+ else _gY->clearAll();
+ if (_gZ == NULL) _gZ = new groupOfElements();
+ else _gZ->clearAll();
+
+ InterpolationOperations::intersectionBetweenTwoFaces(f1234,f1562,*_gX);
+ InterpolationOperations::intersectionBetweenTwoFaces(f1234,f1584,*_gY);
+ InterpolationOperations::intersectionBetweenTwoFaces(f1562,f1584,*_gZ);
this->pbcFEConstraintElementGroup(l34.begin(),l34.end(),_gX,_cVertices,_v4,_v1);
this->pbcFEConstraintElementGroup(l78.begin(),l78.end(),_gX,_cVertices,_v8,_v1);
@@ -1727,29 +1593,6 @@ void pbcConstraintElementGroup::createPeriodicMeshConstraintElementGroup(){
int dim = _solver->getMicroBC()->getDim();
const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- /**fixed node in case of periodic boundary condition and corner nodes does not exist**/
- if ((_cVertices.size() ==0) and (_solver->getMicroBC()->getOrder()==1)){
- if (dim==2){
- MVertex* v = *(bgroup[0]->vbegin());
- _cVertices.insert(v);
- _cVertices.insert(findNearestVertex(bgroup[2],v));
- for (std::set<MVertex*>::iterator it=_cVertices.begin(); it !=_cVertices.end(); it++){
- _bVertices.insert(*it);
- }
- }
- else if (dim==3) {
- MVertex* v = *(bgroup[0]->vbegin());
- _cVertices.insert(v);
- _cVertices.insert(findNearestVertex(bgroup[3],v));
- v = *(bgroup[1]->vbegin());
- _cVertices.insert(v);
- _cVertices.insert(findNearestVertex(bgroup[4],v));
- for (std::set<MVertex*>::iterator it=_cVertices.begin(); it !=_cVertices.end(); it++){
- _bVertices.insert(*it);
- };
- };
- };
-
if (dim ==3){
// periodic of edges
periodicConditionForPeriodicMesh(l56.begin(),l56.end(),l12.begin(),l12.end(),_cVertices,_v5,_v1);
@@ -1787,13 +1630,12 @@ void pbcConstraintElementGroup::createPeriodicMeshConstraintElementGroup(){
}
void pbcConstraintElementGroup::createCubicSplineConstraintElementGroup(){
-
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
int dim = _solver->getMicroBC()->getDim();
-
+ nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
// create basis
- this->createPolynomialBasePoints();
+ std::set<MVertex*> newVirVertices; // new vertices created in interpolation basis
+ InterpolationOperations::createPolynomialBasePoints(pbc,_v1,_v2,_v4,_v5,l12,l41,l15,_xBase,_yBase,_zBase,newVirVertices);
+ _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
if (dim ==2){
this->cubicSplineFormCondition(l12.begin(),l12.end(),_xBase,0);
@@ -1803,6 +1645,9 @@ void pbcConstraintElementGroup::createCubicSplineConstraintElementGroup(){
this->cubicSplinePeriodicCondition(l23.begin(),l23.end(),_yBase,1,_cVertices);
}
else if (dim ==3){
+ // in 3D new ver are b-one
+ _bVertices.insert(newVirVertices.begin(),newVirVertices.end());
+
this->cubicSplineFormCondition(l12.begin(),l12.end(),_xBase,0);
this->cubicSplineFormCondition(l41.begin(),l41.end(),_yBase,1);
this->cubicSplineFormCondition(l15.begin(),l15.end(),_zBase,2);
@@ -1845,7 +1690,9 @@ void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
int dim = _solver->getMicroBC()->getDim();
const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
// create base
- this->createPolynomialBasePoints();
+ std::set<MVertex*> newVirVertices;
+ InterpolationOperations::createPolynomialBasePoints(pbc,_v1,_v2,_v4,_v5,l12,l41,l15,_xBase,_yBase,_zBase,newVirVertices);
+ _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
if (dim ==2){
this->lagrangeFormCondition(l12.begin(),l12.end(),_xBase);
@@ -1893,9 +1740,10 @@ void pbcConstraintElementGroup::createShiftedLagrangeConstraintElementGroup(){
void pbcConstraintElementGroup::createLagrangeConstraintElementGroup(){
int dim = _solver->getMicroBC()->getDim();
nonLinearPeriodicBC* pbc = static_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
// create base
- this->createPolynomialBasePoints();
+ std::set<MVertex*> newVirVertices;
+ InterpolationOperations::createPolynomialBasePoints(pbc,_v1,_v2,_v4,_v5,l12,l41,l15,_xBase,_yBase,_zBase,newVirVertices);
+ _virtualVertices.insert(newVirVertices.begin(),newVirVertices.end());
if (dim ==2){
this->lagrangeFormCondition(l12.begin(),l12.end(),_xBase,_v1,_v1);
@@ -1904,6 +1752,8 @@ void pbcConstraintElementGroup::createLagrangeConstraintElementGroup(){
this->lagrangePeriodicCondition(l23.begin(),l23.end(),_yBase,_cVertices,_v2,_v1);
}
else if (dim ==3){
+ _bVertices.insert(newVirVertices.begin(),newVirVertices.end());
+
this->lagrangeFormCondition(l12.begin(),l12.end(),_xBase,_v1,_v1);
this->lagrangeFormCondition(l41.begin(),l41.end(),_yBase,_v1,_v1);
this->lagrangeFormCondition(l15.begin(),l15.end(),_zBase,_v1,_v1);
@@ -2216,10 +2066,10 @@ void pbcConstraintElementGroup::createMixBCConstraintElementGroupDirectionFollow
constraintElement* ppc = new periodicMeshConstraint(_solver->getMicroBC(), _LagSpace,_LagSpace,_MultSpace,-1,v2);
_allConstraint.insert(ppc);
- // create local basis
- this->getInterpolatioPointsFromExistingVertices(interpDeg[0],v1,v4,bgroup[0]->vbegin(),bgroup[0]->vend(),_xBase);
- this->getInterpolatioPointsFromExistingVertices(interpDeg[1],v1,v2,bgroup[1]->vbegin(),bgroup[1]->vend(),_yBase);
-
+ // 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());
}
@@ -2377,339 +2227,387 @@ void pbcConstraintElementGroup::createCornerConstraintElementGroupShiftedBPC(){
}
void pbcConstraintElementGroup::createCornerConstraintElementGroup(){
- 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
-
+
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());
const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
- if (bgroup.size()>0){
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::CEM or pbc->getPBCMethod() == nonLinearPeriodicBC::PROJECT){
- // for mechanical Dofs
- if (_solver->getMicroBC()->getOrder() == 1){
- // fix all corners mechanical Dofs
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
+ 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");
+ }
+
+ 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);
- }
+ 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());
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC() ,_LagSpace,_MultSpace,mechanics[i],v);
+ _allConstraint.insert(cel);
}
- }
- else if (_solver->getMicroBC()->getOrder() == 2){
- // periodic between
- if (_cVertices.size() >0){
- MVertex* vn = *(_cVertices.begin());
- for(std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- if (it != _cVertices.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);
- }
+ }
+ }
+ 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);
- }
- }
- }
+ 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
- if (_cVertices.size() > 0){
- MVertex* oneV = *(_cVertices.begin());
- for (std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.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);
- }
- }
- }
- }
- 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* 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);
+ }
+ }
+ }
+ }
+ 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
+ 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);
+ _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;
- 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::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],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);
- _allConstraint.insert(ppc);
- }
+ 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 (int i=0; i< compConstitutive.size(); i++){
- constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
- compConstitutive[i], _solver->getDomainVector());
- _allConstraint.insert(celExtra);
- }
+ //
+ // 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);
- }
+ // 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::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 (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);
+ }
+ }
}
- }
- 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;
- for (int i=0; i< mechanics.size(); i++){
- constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],vp1);
- _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;
+ for (int i=0; i< mechanics.size(); i++){
+ constraintElement* cel = new periodicMeshConstraint(_solver->getMicroBC(),_LagSpace,_MultSpace,mechanics[i],vp1);
+ _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 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 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);
+ }
}
-
}
- else if (pbc->getPBCMethod() == nonLinearPeriodicBC::FE_LIN or pbc->getPBCMethod() == nonLinearPeriodicBC::FE_QUA){
+ for (int i=0; i< compConstitutive.size(); i++){
+ constraintElement* celExtra = new supplementConstraintAverageExtraDofValue(_solver->getMicroBC() ,_MultSpace,
+ compConstitutive[i], _solver->getDomainVector());
+ _allConstraint.insert(celExtra);
+ }
- 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);
- }
- }
+ // 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);
+ }
+ }
- 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)){
- activeCornerVertices.insert(vg);
- break;
- }
- }
- }
+ }
+ 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;
+ }
+ }
+ }
- // for mechanics
- if (_solver->getMicroBC()->getOrder() == 1){
- for (std::set<MVertex*>::iterator it = activeCornerVertices.begin(); it!= activeCornerVertices.end(); it++){
+ // 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);
- }
+ 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);
}
+ }
+
+ }
+ 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);
- }
- else if (_solver->getMicroBC()->getOrder() == 2){
- if (_cVertices.size() >0){
- MVertex* vn = *(_cVertices.begin());
- for(std::set<MVertex*>::iterator it = _cVertices.begin(); it!= _cVertices.end(); it++){
- if (it != _cVertices.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() > 2){
- 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 (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 non-constitutive extraDofs
- 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);
- }
+ // 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
- Msg::Error("pbcConstraintElementGroup::createCornerConstraintElementGroup: case undefined");
+
}
+ else
+ Msg::Error("pbcConstraintElementGroup::createCornerConstraintElementGroup: case undefined");
}
};
@@ -2780,8 +2678,9 @@ void pbcConstraintElementGroup::createConstraintGroup(){
if (pbc->isShifted()){
this->createCornerConstraintElementGroupShiftedBPC();
}
- else
+ else{
this->createCornerConstraintElementGroup();
+ }
}
}
else if (_solver->getMicroBC()->getType()==nonLinearMicroBC::LDBC){
@@ -2812,13 +2711,8 @@ void pbcConstraintElementGroup::createConstraintGroup(){
void pbcConstraintElementGroup::clearAllConstraints(){
_allConstraint.clear();
_virtualVertices.clear();
- _pbcNodes.clear();
};
-void pbcConstraintElementGroup::switchMicroBC(){
- this->clearAllConstraints(); // clear all existing constraint elements
- this->createConstraintGroup(); // create new groups
-};
void pbcConstraintElementGroup::createPBCFromFile(FILE* fp, GModel* pModel){
if (fp == NULL) return;
@@ -2876,42 +2770,6 @@ SVector3 pbcConstraintElementGroup::getUniformDisp(MVertex* v){
}
return u;
};
-
-SVector3 pbcConstraintElementGroup::getPerturbation(dofManager<double>* pmanager, MVertex* v){
- SVector3 ubar = this->getUniformDisp(v);
- std::vector<Dof> keys;
- std::vector<double> val;
- std::vector<int> constrainedComp(_solver->getMicroBC()->getConstrainedComps().begin(), _solver->getMicroBC()->getConstrainedComps().end());
- getKeysFromVertex(_LagSpace,v,constrainedComp,keys);
- pmanager->getDofValue(keys,val);
- ubar *= -1.;
- for (int i=0; i<keys.size(); i++){
- ubar[i] += val[i];
- }
- return ubar;
-};
-
-SVector3 pbcConstraintElementGroup::getTangentCubicSpline(dofManager<double>* pmanager, MVertex* v, const int dir){
- SVector3 vec(0.,0.,0.);
- std::vector<int> constrainedComp(_solver->getMicroBC()->getConstrainedComps().begin(), _solver->getMicroBC()->getConstrainedComps().end());
- if (_solver->getMicroBC()->getType() != nonLinearMicroBC::PBC) {
- Msg::Warning("this used for periodic BC and cubic spline method");
- }
- else {
- nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
- if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM){
- std::vector<Dof> dofs;
- std::vector<double> vals;
- cubicSplineConstraintElement::addKeysToVertex(_LagSpace,v,constrainedComp,dofs,dir);
- pmanager->getDofValue(dofs,vals);
- for (int i=0; i<vals.size(); i++){
- vec[i] += vals[i];
- }
- }
- }
- return vec;
-};
-
void pbcConstraintElementGroup::cubicSplineToFile(std::vector<MVertex*>& vlist, dofManager<double>* pmanager, FILE* file, const int dir){
if (vlist.size() ==0) return;
int N= 100;
@@ -2985,6 +2843,43 @@ void pbcConstraintElementGroup::lagrangeToFile(std::vector<MVertex*>& vlist, dof
}
};
+SVector3 pbcConstraintElementGroup::getPerturbation(dofManager<double>* pmanager, MVertex* v){
+ SVector3 ubar = this->getUniformDisp(v);
+ std::vector<Dof> keys;
+ std::vector<double> val;
+ std::vector<int> constrainedComp(_solver->getMicroBC()->getConstrainedComps().begin(), _solver->getMicroBC()->getConstrainedComps().end());
+ getKeysFromVertex(_LagSpace,v,constrainedComp,keys);
+ pmanager->getDofValue(keys,val);
+ ubar *= -1.;
+ for (int i=0; i<keys.size(); i++){
+ ubar[i] += val[i];
+ }
+ return ubar;
+};
+
+SVector3 pbcConstraintElementGroup::getTangentCubicSpline(dofManager<double>* pmanager, MVertex* v, const int dir){
+ SVector3 vec(0.,0.,0.);
+ std::vector<int> constrainedComp(_solver->getMicroBC()->getConstrainedComps().begin(), _solver->getMicroBC()->getConstrainedComps().end());
+ if (_solver->getMicroBC()->getType() != nonLinearMicroBC::PBC) {
+ Msg::Warning("this used for periodic BC and cubic spline method");
+ }
+ else {
+ nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
+ if (pbc->getPBCMethod() == nonLinearPeriodicBC::CSIM){
+ std::vector<Dof> dofs;
+ std::vector<double> vals;
+ cubicSplineConstraintElement::addKeysToVertex(_LagSpace,v,constrainedComp,dofs,dir);
+ pmanager->getDofValue(dofs,vals);
+ for (int i=0; i<vals.size(); i++){
+ vec[i] += vals[i];
+ }
+ }
+ }
+ return vec;
+};
+
+
+
void pbcConstraintElementGroup::writePertBoundaryToFile(dofManager<double>* pmanager, const int e, const int g, const int iter ){
if (_solver->getMicroBC()->getType() == nonLinearMicroBC::PBC){
nonLinearPeriodicBC* pbc = dynamic_cast<nonLinearPeriodicBC*>(_solver->getMicroBC());
@@ -3167,17 +3062,6 @@ void pbcConstraintElementGroup::applyLinearConstraintsToSystem(dofManager<double
};
};
-void pbcConstraintElementGroup::setLocalBasis(const SVector3& n, const SVector3& t, const SVector3& b){
- _e1 = n; _e1.normalize();
- _e2 = t; _e2.normalize();
- _e3 = b; _e3.normalize();
-};
-
-void pbcConstraintElementGroup::setRotatedBasis(const SVector3& n, const SVector3& t, const SVector3 b){
- _e1Rotated = n; _e1Rotated.normalize();
- _e2Rotated = t; _e2Rotated.normalize();
- _e3Rotated = b; _e3Rotated.normalize();
-};
SPoint3 pbcConstraintElementGroup::getRootPoint(){
if (_v1 != NULL) return _v1->point();
@@ -3185,7 +3069,7 @@ SPoint3 pbcConstraintElementGroup::getRootPoint(){
const std::vector<groupOfElements*>& bgroup = _solver->getMicroBC()->getBoundaryGroupOfElements();
if (_solver->getMicroBC()->getDim() == 2){
std::set<MVertex*> vv;
- groupIntersection(bgroup[0],bgroup[1],vv);
+ InterpolationOperations::groupIntersection(bgroup[0],bgroup[1],vv);
if (vv.size() == 1){
_v1 = *(vv.begin());
}
diff --git a/NonLinearSolver/periodicBC/pbcCreateConstraints.h b/NonLinearSolver/periodicBC/pbcCreateConstraints.h
index 273f696486d634cec307952db62d5c1705cf86bd..828895456d82e48bc8b2302adcc0ce13fa9ba1aa 100644
--- a/NonLinearSolver/periodicBC/pbcCreateConstraints.h
+++ b/NonLinearSolver/periodicBC/pbcCreateConstraints.h
@@ -72,6 +72,47 @@ class GModel;
1--------l12---------2
*/
+namespace InterpolationOperations{
+ void groupIntersection(groupOfElements* g1, groupOfElements* g2, std::set<MVertex*>& ver);
+ void groupIntersection(std::set<MVertex*>& g1, std::set<MVertex*>& g2, std::set<MVertex*>& ver);
+ /* create group of line at intersection of two faces*/
+ void intersectionBetweenTwoFaces(const groupOfElements* face1, const groupOfElements* face2, groupOfElements& line);
+ /* get interpolation points on a line*/
+ void createLinePolynomialBase(const int degree, // degree
+ const bool createNewVer, // true if new vertices are created false if interpolation base is taken from existing vertices
+ MVertex* vstart, MVertex* vend, // vstart and vend
+ std::set<MVertex*>& line, // current line
+ std::vector<MVertex*>& base, // base obtained
+ std::set<MVertex*>& virtualVertices // virtualVertices to store all virtual virtices which has been created
+ );
+ /*
+ This function allows to find a nearest vertex in a group vertex to one vertex
+ by using the minimal distance flag
+ */
+ template<class Iterator>
+ MVertex* findNearestVertex(Iterator itbegin, Iterator itend, MVertex* vp);
+ MVertex* findNearestVertex(groupOfElements* g, MVertex* vp);
+
+
+ void createGroupLinearLineElementFromVerticesVector(std::vector<MVertex*>& b, groupOfElements& gr);
+ void createGroupQuadraticLineElementFromVerticesVector(std::vector<MVertex*>& b, groupOfElements& gr, std::set<MVertex*>& virVertices);
+
+ /*create a surface group to interpolate form two boundaries for FE method */
+ void createLinearElementBCGroup3DFromTwoEdges(std::vector<MVertex*>& vlx, std::vector<MVertex*>& vly,
+ MVertex* vc, groupOfElements& g, std::set<MVertex*>& newVer);
+ void createQuadraticElementBCGroup3DFromTwoEdges(std::vector<MVertex*>& vlx, std::vector<MVertex*>& vly,
+ MVertex* vc, groupOfElements& g, std::set<MVertex*>& newVer);
+
+ void createPolynomialBasePoints(nonLinearPeriodicBC* pbc,
+ MVertex* v1, MVertex* v2, MVertex* v4, MVertex* v5,
+ std::set<MVertex*>& l12, std::set<MVertex*>& l14,std::set<MVertex*>& l15,
+ std::vector<MVertex*>& xbase, std::vector<MVertex*>& ybase, std::vector<MVertex*>& zbase,
+ std::set<MVertex*>& newVer);
+ bool isInside(MVertex* v, MElement* e);
+
+};
+
+
class groupOfElementsDomainDecomposition{
public:
int physical;
@@ -90,59 +131,36 @@ class pbcConstraintElementGroup{
FunctionSpaceBase* _LagSpace, *_MultSpace; // lagspace and mult space
std::set<constraintElement*> _allConstraint; // all constraint
+
std::set<MVertex*> _cVertices; // corner vertices
std::set<MVertex*> _bVertices; // edge vertices, in 2D corner and edger vertices are same
std::set<MVertex*> _virtualVertices; // new vertices added to systems for polynomial interpolation
- std::vector<MVertex*> _xBase, _yBase, _zBase; // interpolation basis
+
groupOfElements* _g;
/*in case of 2D*/
- groupOfElements* f1234; // surface
- std::set<MVertex*> l12, l23, l34,l41; // 4 edges
+ std::set<MVertex*> l12, l23, l34,l41; // 4 edges
//add surfaces of RVE in 3D cases
- groupOfElements *f5678, *f1562, *f4873, *f1584, *f2673;
+ groupOfElements * f1234, *f5678, *f1562, *f4873, *f1584, *f2673;
// add 8 edges of RVE in 3D case
std::set<MVertex*> l56, l67, l78, l85, l15, l26, l37, l48;
MVertex* _v1, *_v2, *_v4, *_v5; // all control points
MVertex* _v3, *_v6, *_v7, *_v8;
- groupOfElements* _gX, *_gY, *_gZ; // interpolation edges
- groupOfElements* _gXY, *_gYZ, *_gZX; // interpolation faces
-
- SVector3 _periodic12, _periodic14, _periodic15; //
- SVector3 _e1, _e2, _e3; // local basis
- SVector3 _e1Rotated, _e2Rotated, _e3Rotated; // rotate basis
-
- // 6 periodic nodes in case of missing corner nodes
- std::vector<MVertex*> _pbcNodes;
+
+ SVector3 _periodic12, _periodic14, _periodic15; // periodicity computed from RVE
+
+ std::vector<MVertex*> _xBase, _yBase, _zBase; // default interpolation basis
+ groupOfElements* _gX, *_gY, *_gZ, *_gXY, *_gYZ, *_gZX; // basis PC based on FE method
std::vector<groupOfElementsDomainDecomposition*> _boundaryDGMap;
std::map<partDomain*,FunctionSpaceBase*> _allDGMultSpace; // all multspace for DG case
- protected:
- void groupIntersection(groupOfElements* g1, groupOfElements* g2, std::set<MVertex*>& ver) const;
- void groupIntersection(groupOfElements* g1, groupOfElements* g2, groupOfElements* &g) const;
- void groupIntersection(std::set<MVertex*>& g1, std::set<MVertex*>& g2, std::set<MVertex*>& ver) const;
- /* create group of line at intersection of two faces*/
- void intersectionBetweenTwoFaces(const groupOfElements* face1, const groupOfElements* face2, groupOfElements* &line) const;
- void intersectionBetweenTwoLines(const groupOfElements* line1, const groupOfElements* line2, groupOfElements* & pt) const;
+ private:
+
+
void __init__();
- // add control point in case of missing
- void addControlPoints();
- // add other corner nodes
- void addOtherCornerPoints();
- // create interpolation basis
- void createPolynomialBasePoints();
- /*create a surface group to interpolate form two boundaries for FE method */
- void createBCGroup3D(std::vector<MVertex*>& vlx, std::vector<MVertex*>& vly,
- MVertex* vc, groupOfElements* g);
- /*create all interpolate surfaces for FE method only*/
- void createAllBCGroups();
- /*
- This function allows to find a nearest vertex in a group vertex to one vertex
- by using the minimal distance flag
- */
- template<class Iterator>
- MVertex* findNearestVertex(Iterator itbegin, Iterator itend, MVertex* vp);
- MVertex* findNearestVertex(groupOfElements* g, MVertex* vp);
+
+ // in this case, vinit and vend using for constructing the distance function
+
/*
Enforce periodic boundary condition on conformal meshes
positive part: posivtiveitbegin-->positiveitend
@@ -159,15 +177,7 @@ class pbcConstraintElementGroup{
Iterator negativeitbegin, Iterator negativeitend,
std::set<MVertex*>& others, MVertex* vrootP=NULL, MVertex* vrootN=NULL);
- /* get interpolation points on a line*/
- void createNewVerticesForInterpolationPoints(int degree, MVertex* vinit, MVertex* vend,
- std::vector<MVertex*>& base,
- std::set<MVertex*>& others);
- // in this case, vinit and vend using for constructing the distance function
- template<class Iterator>
- void getInterpolatioPointsFromExistingVertices(int degree, MVertex* vinit, MVertex* vend,
- Iterator itbegin, Iterator itend,
- std::vector<MVertex*>& base);
+
// for cubic spline
template<class Iterator>
@@ -216,8 +226,11 @@ class pbcConstraintElementGroup{
void lagrangePeriodicConditionCoonsPatch(Iterator itbegin, Iterator itend, MVertex* vref,
std::vector<MVertex*>& xlist,std::vector<MVertex*>& ylist,
std::set<MVertex*>& others ,MVertex* vrootP=NULL, MVertex* vrootN = NULL);
-
- bool isInside(MVertex* v, MElement* e) const;
+
+ template<class Iterator>
+ void pbcFormFEConstraintElementGroup(Iterator itbegin, Iterator itend, groupOfElements* g,
+ std::set<MVertex*>& others, MVertex* vrootP=NULL, MVertex* vrootN=NULL);
+
template<class Iterator>
void pbcFEConstraintElementGroup(Iterator itbegin, Iterator itend, groupOfElements* g,
std::set<MVertex*>& others, MVertex* vrootP=NULL, MVertex* vrootN=NULL);
@@ -248,6 +261,7 @@ class pbcConstraintElementGroup{
SVector3 getUniformDisp(MVertex* v);
SVector3 getPerturbation(dofManager<double>* pmanager, MVertex* v);
SVector3 getTangentCubicSpline(dofManager<double>* pmanager, MVertex* v,const int dir);
+ // write perturbation to FILE
void cubicSplineToFile(std::vector<MVertex*>& vlist, dofManager<double>* pmanager, FILE* file, const int dir);
void lagrangeToFile(std::vector<MVertex*>& vlist, dofManager<double>* pmanager, FILE* file);
@@ -256,7 +270,6 @@ class pbcConstraintElementGroup{
FunctionSpaceBase* lagspace, FunctionSpaceBase* multspace);
virtual ~pbcConstraintElementGroup();
void createConstraintGroup();
- void switchMicroBC();
void clearAllConstraints();
std::set<constraintElement*>::iterator constraintGroupBegin(){return _allConstraint.begin();};
std::set<constraintElement*>::iterator constraintGroupEnd(){return _allConstraint.end();};
@@ -280,8 +293,6 @@ class pbcConstraintElementGroup{
void numberDof_virtualVertices(dofManager<double>* pmanager);
void applyLinearConstraintsToSystem(dofManager<double>* pmanager);
- void setLocalBasis(const SVector3& n, const SVector3& t, const SVector3& b);
- void setRotatedBasis(const SVector3& n, const SVector3& t, const SVector3 b);
SPoint3 getRootPoint();
};