diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 4f68cb58a30a4bf529421d97557e5fdf70c530d1..e1577642b6fd499afe2a5c029723e14997d2a15c 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -148,28 +148,37 @@ dgGroupOfElements::~dgGroupOfElements(){
delete _elementVolume;
}
+void dgGroupOfConnections::addElement(int iElement, int iClosure) {
+ _elementId.push_back(iElement);
+ _closuresId.push_back(iClosure);
+}
+
+
+
void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
- // compute all closures
- // compute closures for the interpolation
- _left.push_back(iElLeft);
- MElement &elLeft = *_groupLeft.getElement(iElLeft);
- int ithFace, sign, rot;
- elLeft.getFaceInfo (topoFace, ithFace, sign, rot);
- int closureIdLeft = _fsLeft->getFaceClosureId(ithFace, sign, rot);
- _closuresIdLeft.push_back(closureIdLeft);
- if(iElRight>=0){
- _right.push_back(iElRight);
- MElement &elRight = *_groupRight.getElement(iElRight);
- elRight.getFaceInfo (topoFace, ithFace, sign, rot);
- _closuresIdRight.push_back(_fsRight->getFaceClosureId(ithFace, sign, rot));
- }
// compute the face element that correspond to the geometrical closure
// get the vertices of the face
std::vector<MVertex*> vertices;
- const std::vector<int> & geomClosure = elLeft.getFunctionSpace()->getFaceClosure(closureIdLeft);
+ int closureId;
+ const MElement *el;
+ int ithFace, sign, rot;
+ if(iElRight>=0){
+ el = _connections[1]->getGroupOfElements().getElement(iElRight);
+ el->getFaceInfo(topoFace, ithFace, sign, rot);
+ closureId = _connections[1]->getFunctionSpace()->getFaceClosureId(ithFace,sign,rot);
+ _connections[1]->addElement(iElRight, closureId);
+ }
+
+ el = _connections[0]->getGroupOfElements().getElement(iElLeft);
+ el->getFaceInfo(topoFace, ithFace, sign, rot);
+ closureId = _connections[0]->getFunctionSpace()->getFaceClosureId(ithFace,sign,rot);
+ _connections[0]->addElement(iElLeft, closureId);
+
+ const std::vector<int> & geomClosure = el->getFunctionSpace()->getFaceClosure(closureId);
for (int j=0; j<geomClosure.size() ; j++)
- vertices.push_back( elLeft.getVertex(geomClosure[j]) );
+ vertices.push_back( const_cast<MElement*>(el)->getVertex(geomClosure[j]) );
+
// triangular face
if (topoFace.getNumVertices() == 3){
switch(vertices.size()){
@@ -180,8 +189,7 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
case 21 : _faces.push_back(new MTriangleN (vertices,5) ); break;
default : throw;
}
- }
- else if (topoFace.getNumVertices() == 4){
+ } else if (topoFace.getNumVertices() == 4){
switch(vertices.size()){
case 4 : _faces.push_back(new MQuadrangle (vertices) ); break;
case 8 : _faces.push_back(new MQuadrangle8 (vertices) ); break;
@@ -191,50 +199,50 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
default : throw;
}
}
- // quad face 2 do
else {
throw;
}
}
void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight) {
- _left.push_back(iElLeft);
- MElement &elLeft = *_groupLeft.getElement(iElLeft);
- int ithEdge, sign;
- elLeft.getEdgeInfo (topoEdge, ithEdge, sign);
- _closuresIdLeft.push_back(_fsLeft->getEdgeClosureId(ithEdge, sign));
+ std::vector<MVertex*> vertices;
+ const MElement *el ;
+ int closureId;
+ int ithFace, sign;
if(iElRight>=0) {
- _right.push_back(iElRight);
- _groupRight.getElement(iElRight)->getEdgeInfo (topoEdge, ithEdge, sign);
- _closuresIdRight.push_back(_fsRight->getEdgeClosureId(ithEdge,sign));
+ el = _connections[1]->getGroupOfElements().getElement(iElRight);
+ el->getEdgeInfo(topoEdge, ithFace, sign);
+ closureId = _connections[1]->getFunctionSpace()->getEdgeClosureId(ithFace,sign);
+ _connections[1]->addElement(iElRight, closureId);
}
- std::vector<MVertex*> vertices;
- const std::vector<int> &geomClosure = elLeft.getFunctionSpace()->getEdgeClosure(_closuresIdLeft.back());
+
+ el = _connections[0]->getGroupOfElements().getElement(iElLeft);
+ el->getEdgeInfo(topoEdge, ithFace, sign);
+ closureId = _connections[0]->getFunctionSpace()->getEdgeClosureId(ithFace,sign);
+ _connections[0]->addElement(iElLeft, closureId);
+ const std::vector<int> & geomClosure = el->getFunctionSpace()->getEdgeClosure(closureId);
for (int j=0; j<geomClosure.size() ; j++)
- vertices.push_back( elLeft.getVertex(geomClosure[j]) );
+ vertices.push_back( const_cast<MElement*>(el)->getVertex(geomClosure[j]) );
switch(vertices.size())
- {
+ {
case 2 : _faces.push_back(new MLine (vertices) ); break;
case 3 : _faces.push_back(new MLine3 (vertices) ); break;
default : _faces.push_back(new MLineN (vertices) ); break;
- }
+ }
}
void dgGroupOfFaces::addVertex(MVertex *topoVertex, int iElLeft, int iElRight){
- _left.push_back(iElLeft);
- MElement &elLeft = *_groupLeft.getElement(iElLeft);
int ithVertex;
- elLeft.getVertexInfo (topoVertex, ithVertex);
- _closuresIdLeft.push_back(ithVertex);
+ _connections[0]->getGroupOfElements().getElement(iElLeft)->getVertexInfo(topoVertex, ithVertex);
+ _connections[0]->addElement(iElLeft,ithVertex);
if(iElRight>=0){
- _right.push_back(iElRight);
- MElement &elRight = *_groupRight.getElement(iElRight);
- elRight.getVertexInfo (topoVertex, ithVertex);
- _closuresIdRight.push_back(ithVertex);
+ _connections[1]->getGroupOfElements().getElement(iElRight)->getVertexInfo(topoVertex, ithVertex);
+ _connections[1]->addElement(iElRight, ithVertex);
}
_faces.push_back(new MPoint(topoVertex) );
}
+
void dgGroupOfFaces::init(int pOrder) {
if (!_faces.size())return;
_fsFace = _faces[0]->getFunctionSpace (pOrder);
@@ -243,12 +251,10 @@ void dgGroupOfFaces::init(int pOrder) {
_collocation = new fullMatrix<double> (_integration->size1(), _fsFace->coefficients.size1());
_detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
_interfaceSurface = new fullMatrix<double>(_faces.size(),1);
- for (size_t i=0; i<_closuresLeft.size(); i++)
- _integrationPointsLeft.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,_closuresLeft[i]));
- for (size_t i=0; i<_closuresRight.size(); i++)
- _integrationPointsRight.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsRight,_closuresRight[i]));
+ for(int i=0; i<_connections.size(); i++) {
+ _connections[i]->init();
+ }
double f[256];
-
for (int j=0;j<_integration->size1();j++) {
_fsFace->f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
const double weight = (*_integration)(j,3);
@@ -257,7 +263,6 @@ void dgGroupOfFaces::init(int pOrder) {
(*_collocation)(j,k) = f[k];
}
}
-
for (int i=0;i<_faces.size();i++){
MElement *f = _faces[i];
for (int j=0;j< _integration->size1() ; j++ ){
@@ -266,70 +271,6 @@ void dgGroupOfFaces::init(int pOrder) {
(*_interfaceSurface)(i,0) += (*_integration)(j,3)*(*_detJac)(j,i);
}
}
-
- // compute data on quadrature points : normals and dPsidX
- double g[256][3];
- _normals = new fullMatrix<double> (3,_integration->size1()*_faces.size());
- _dPsiLeftDxOnQP = new fullMatrix<double> ( _integration->size1()*3,_fsLeft->points.size1()*_faces.size());
- if(_fsRight){
- _dPsiRightDxOnQP = new fullMatrix<double> ( _integration->size1()*3,_fsRight->points.size1()*_faces.size());
- }
- int index = 0;
- for (size_t i=0; i<_faces.size();i++){
-
- const std::vector<int> &closureLeft=getClosureLeft(i);
- fullMatrix<double> &intLeft=_integrationPointsLeft[_closuresIdLeft[i]];
- double jac[3][3],ijac[3][3];
- for (int j=0; j<intLeft.size1(); j++){
- _fsLeft->df((intLeft)(j,0),(intLeft)(j,1),(intLeft)(j,2),g);
- getElementLeft(i)->getJacobian ((intLeft)(j,0), (intLeft)(j,1), (intLeft)(j,2), jac);
- inv3x3(jac,ijac);
- //compute dPsiLeftDxOnQP
- // (iQP*3+iXYZ , iFace*NPsi+iPsi)
- int nPsi = _fsLeft->coefficients.size1();
- for (int iPsi=0; iPsi< nPsi; iPsi++) {
- (*_dPsiLeftDxOnQP)(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
- (*_dPsiLeftDxOnQP)(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
- (*_dPsiLeftDxOnQP)(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
- }
- //compute face normals
- double &nx=(*_normals)(0,index);
- double &ny=(*_normals)(1,index);
- double &nz=(*_normals)(2,index);
- double nu=0,nv=0,nw=0;
- for (size_t k=0; k<closureLeft.size(); k++){
- int inode=closureLeft[k];
- nu += g[inode][0];
- nv += g[inode][1];
- nw += g[inode][2];
- }
- nx = nu*ijac[0][0]+nv*ijac[0][1]+nw*ijac[0][2];
- ny = nu*ijac[1][0]+nv*ijac[1][1]+nw*ijac[1][2];
- nz = nu*ijac[2][0]+nv*ijac[2][1]+nw*ijac[2][2];
- double norm = sqrt(nx*nx+ny*ny+nz*nz);
- nx/=norm;
- ny/=norm;
- nz/=norm;
- index++;
- }
- // there is nothing on the right for boundary groups
- if(_fsRight){
- fullMatrix<double> &intRight=_integrationPointsRight[_closuresIdRight[i]];
- for (int j=0; j<intRight.size1(); j++){
- _fsRight->df((intRight)(j,0),(intRight)(j,1),(intRight)(j,2),g);
- getElementRight(i)->getJacobian ((intRight)(j,0), (intRight)(j,1), (intRight)(j,2), jac);
- inv3x3(jac,ijac);
- //compute dPsiRightDxOnQP
- // (iQP*3+iXYZ , iFace*NPsi+iPsi)
- int nPsi = _fsRight->coefficients.size1();
- for (int iPsi=0; iPsi< nPsi; iPsi++) {
- (*_dPsiRightDxOnQP)(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
- (*_dPsiRightDxOnQP)(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
- (*_dPsiRightDxOnQP)(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
- }
- }
- }
- }
}
dgGroupOfFaces::~dgGroupOfFaces()
@@ -338,23 +279,17 @@ dgGroupOfFaces::~dgGroupOfFaces()
delete _redistribution;
delete _collocation;
delete _detJac;
- delete _normals;
- delete _dPsiLeftDxOnQP;
- delete _dPsiRightDxOnQP;
delete _interfaceSurface;
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices):
- _groupLeft(elGroup),_groupRight(elGroup)
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices)
{
+ _connections.push_back(new dgGroupOfConnections(elGroup, *this, pOrder));
_boundaryTag=boundaryTag;
if(boundaryTag==""){
Msg::Warning ("empty boundary tag, group of boundary faces not created");
return;
}
- _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
- _closuresLeft = _fsLeft->vertexClosure;
- _fsRight=NULL;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumVertices(); j++){
@@ -366,17 +301,14 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
init(pOrder);
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges):
- _groupLeft(elGroup),_groupRight(elGroup)
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges)
{
+ _connections.push_back(new dgGroupOfConnections(elGroup,*this, pOrder));
_boundaryTag=boundaryTag;
if(boundaryTag==""){
Msg::Warning ("empty boundary tag, group of boundary faces not created");
return;
}
- _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
- _closuresLeft = _fsLeft->edgeClosure;
- _fsRight=NULL;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumEdges(); j++){
@@ -388,15 +320,12 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
init(pOrder);
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces):
- _groupLeft(elGroup),_groupRight(elGroup)
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces)
{
+ _connections.push_back(new dgGroupOfConnections(elGroup, *this, pOrder));
_boundaryTag=boundaryTag;
if(boundaryTag=="")
throw;
- _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
- _closuresLeft = _fsLeft->faceClosure;
- _fsRight=NULL;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumFaces(); j++){
@@ -408,16 +337,20 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
init(pOrder);
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, int numVertices):
- _groupLeft(elGroup),_groupRight(elGroup)
+dgGroupOfConnections::dgGroupOfConnections(const dgGroupOfElements &elementGroup, const dgGroupOfFaces &faceGroup, int pOrder) :
+ _elementGroup(elementGroup),
+ _faceGroup(faceGroup)
+{
+ _fs = _elementGroup.getElement(0)->getFunctionSpace(pOrder);
+}
+
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, int numVertices)
{
- _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
- _fsRight=_groupRight.getElement(0)->getFunctionSpace(pOrder);
- switch (_groupLeft.getElement(0)->getDim()) {
+ _connections.push_back(new dgGroupOfConnections(elGroup, *this, pOrder));
+ _connections.push_back(new dgGroupOfConnections(elGroup, *this, pOrder));
+ switch (_connections[0]->getGroupOfElements().getElement(0)->getDim()) {
case 1 : {
std::map<MVertex*,int> vertexMap;
- _closuresLeft = _fsLeft->vertexClosure;
- _closuresRight = _fsRight->vertexClosure;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumVertices(); j++){
@@ -433,8 +366,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, in
}
case 2 : {
std::map<MEdge,int,Less_Edge> edgeMap;
- _closuresLeft = _fsLeft->edgeClosure;
- _closuresRight = _fsRight->edgeClosure;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumEdges(); j++){
@@ -450,8 +381,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, in
}
case 3 : {
std::map<MFace,int,Less_Face> faceMap;
- _closuresLeft = _fsLeft->faceClosure;
- _closuresRight = _fsRight->faceClosure;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumFaces(); j++){
@@ -472,16 +401,13 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder, in
init(pOrder);
}
-dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder, int numVertices):
- _groupLeft(elGroup1),_groupRight(elGroup2)
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder, int numVertices)
{
- _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
- _fsRight=_groupRight.getElement(0)->getFunctionSpace(pOrder);
- switch (_groupLeft.getElement(0)->getDim()) {
+ _connections.push_back(new dgGroupOfConnections(elGroup1, *this, pOrder));
+ _connections.push_back(new dgGroupOfConnections(elGroup2, *this, pOrder));
+ switch (_connections[0]->getGroupOfElements().getElement(0)->getDim()) {
case 1 : {
std::map<MVertex*,int> vertexMap1;
- _closuresLeft = _fsLeft->vertexClosure;
- _closuresRight = _fsRight->vertexClosure;
for(int i=0; i<elGroup1.getNbElements(); i++){
MElement &el = *elGroup1.getElement(i);
for (int j=0; j<el.getNumVertices(); j++){
@@ -508,8 +434,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
break;
case 2 : {
std::map<MEdge,int,Less_Edge> edgeMap;
- _closuresLeft = _fsLeft->edgeClosure;
- _closuresRight = _fsRight->edgeClosure;
for(int i=0; i<elGroup1.getNbElements(); i++){
MElement &el = *elGroup1.getElement(i);
for (int j=0; j<el.getNumEdges(); j++){
@@ -527,12 +451,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
MEdge edge = el.getEdge(j);
std::map<MEdge,int,Less_Edge>::iterator it = edgeMap.find(edge);
if(it != edgeMap.end()){
-// MElement &el2 = *elGroup1.getElement(it->second);
-// printf("%d %d \n",edge.getVertex(0)->getNum(),edge.getVertex(1)->getNum());
-// for (int k=0;k<el.getNumVertices();k++)printf("%d ",el.getVertex(k)->getNum());
-// printf("\n");
-// for (int k=0;k<el2.getNumVertices();k++)printf("%d ",el2.getVertex(k)->getNum());
-// printf("\n");
addEdge(edge,it->second,i);
}
}
@@ -542,8 +460,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
}
case 3 : {
std::map<MFace,int,Less_Face> faceMap;
- _closuresLeft = _fsLeft->faceClosure;
- _closuresRight = _fsRight->faceClosure;
for(int i=0; i<elGroup1.getNbElements(); i++){
MElement &el = *elGroup1.getElement(i);
for (int j=0; j<el.getNumFaces(); j++){
@@ -583,7 +499,7 @@ void dgGroupOfFaces::mapToInterface ( int nFields,
const std::vector<int> &closureLeft = getClosureLeft(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++){
- v(j, i*nFields + iField) = vLeft(closureLeft[j], _left[i]*nFields + iField);
+ v(j, i*nFields + iField) = vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField);
}
}
}
@@ -594,10 +510,10 @@ void dgGroupOfFaces::mapToInterface ( int nFields,
const std::vector<int> &closureLeft = getClosureLeft(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++){
- v(j, i*2*nFields + iField) = vLeft(closureLeft[j], _left[i]*nFields + iField);
+ v(j, i*2*nFields + iField) = vLeft(closureLeft[j],_connections[0]->getElementId(i)*nFields + iField);
}
for(size_t j =0; j < closureRight.size(); j++)
- v(j, (i*2+1)*nFields + iField) = vRight(closureRight[j], _right[i]*nFields + iField);
+ v(j, (i*2+1)*nFields + iField) = vRight(closureRight[j],_connections[1]->getElementId(i)*nFields + iField);
}
}
}
@@ -614,7 +530,7 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
const std::vector<int> &closureLeft = getClosureLeft(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++)
- vLeft(closureLeft[j], _left[i]*nFields + iField) += v(j, i*nFields + iField);
+ vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*nFields + iField);
}
}
}else{
@@ -623,9 +539,9 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
const std::vector<int> &closureLeft = getClosureLeft(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++)
- vLeft(closureLeft[j], _left[i]*nFields + iField) += v(j, i*2*nFields + iField);
+ vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*2*nFields + iField);
for(size_t j =0; j < closureRight.size(); j++)
- vRight(closureRight[j], _right[i]*nFields + iField) += v(j, (i*2+1)*nFields + iField);
+ vRight(closureRight[j], _connections[1]->getElementId(i)*nFields + iField) += v(j, (i*2+1)*nFields + iField);
}
}
}
@@ -640,7 +556,7 @@ void dgGroupOfFaces::mapLeftFromInterface ( int nFields,
const std::vector<int> &closureLeft = getClosureLeft(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureLeft.size(); j++)
- vLeft(closureLeft[j], _left[i]*nFields + iField) += v(j, i*2*nFields + iField);
+ vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*2*nFields + iField);
}
}
}
@@ -654,31 +570,69 @@ void dgGroupOfFaces::mapRightFromInterface ( int nFields,
const std::vector<int> &closureLeft = getClosureLeft(i);
for (int iField=0; iField<nFields; iField++){
for(size_t j =0; j < closureRight.size(); j++)
- vRight(closureRight[j], _right[i]*nFields + iField) += v(j, (i*2+1)*nFields + iField);
+ vRight(closureRight[j], _connections[1]->getElementId(i)*nFields + iField) += v(j, (i*2+1)*nFields + iField);
}
}
}
-/*
-void dgAlgorithm::buildGroupsOfFaces(GModel *model, int dim, int order,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups){
-}
-
-void dgAlgorithm::buildMandatoryGroups(dgGroupOfElements &eGroup,
- std::vector<dgGroupOfElements*> &partitionedGroups){
-}
+void dgGroupOfConnections::init() {
+ switch(getElement(0)->getDim()) {
+ case 1 : _closures = _fs->vertexClosure; break;
+ case 2 : _closures = _fs->edgeClosure; break;
+ case 3 : _closures = _fs->faceClosure; break;
+ }
+ for (size_t i=0; i<_closures.size(); i++)
+ _integrationPoints.push_back(dgGetFaceIntegrationRuleOnElement(
+ _faceGroup.getPolynomialBasis(), _faceGroup.getIntegrationPointsMatrix(),
+ _fs,_closures[i]));
-void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
- std::vector<dgGroupOfElements*> &partitionedGroups){
+ int nPsi = _fs->points.size1();
+ int nQP = _integrationPoints[0].size1();
+ int size = _elementId.size();
+ // compute data on quadrature points : normals and dPsidX
+ double g[256][3];
+ _normals = fullMatrix<double> (3, nQP*size);
+ _dPsiDxOnQP = fullMatrix<double> ( nQP*3, nPsi*size);
+ int index = 0;
+ for (size_t i=0; i<size;i++){
+ const std::vector<int> &closure = getClosure(i);
+ const fullMatrix<double> &integration = getIntegrationPointsOnElement(i);
+ double jac[3][3],ijac[3][3];
+ for (int j=0; j<integration.size1(); j++){
+ _fs->df(integration(j,0), integration(j,1), integration(j,2),g);
+ getElement(i)->getJacobian (integration(j,0), integration(j,1), integration(j,2), jac);
+ inv3x3(jac,ijac);
+ //compute dPsiDxOnQP
+ // (iQP*3+iXYZ , iFace*NPsi+iPsi)
+ for (int iPsi=0; iPsi< nPsi; iPsi++) {
+ _dPsiDxOnQP(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
+ _dPsiDxOnQP(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
+ _dPsiDxOnQP(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
+ }
+ //compute face normals
+ double &nx=_normals(0,index);
+ double &ny=_normals(1,index);
+ double &nz=_normals(2,index);
+ double nu=0,nv=0,nw=0;
+ for (size_t k=0; k<closure.size(); k++){
+ int inode=closure[k];
+ nu += g[inode][0];
+ nv += g[inode][1];
+ nw += g[inode][2];
+ }
+ nx = nu*ijac[0][0]+nv*ijac[0][1]+nw*ijac[0][2];
+ ny = nu*ijac[1][0]+nv*ijac[1][1]+nw*ijac[1][2];
+ nz = nu*ijac[2][0]+nv*ijac[2][1]+nw*ijac[2][2];
+ double norm = sqrt(nx*nx+ny*ny+nz*nz);
+ nx/=norm;
+ ny/=norm;
+ nz/=norm;
+ index++;
+ }
+ }
}
-*/
-
-//static void partitionGroup (std::vector<MElement *>,){
-//}
// this function creates groups of elements
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 042a48c909e776d321bd5ea974ca498b819984b4..1d65932a958d295db8dc00bde7aff3700e6ccc40 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -123,19 +123,45 @@ public:
friend class dgGroupCollection;
};
+class dgGroupOfFaces;
+
+class dgGroupOfConnections {
+ std::vector<std::vector<int> > _closures;
+ std::vector<int> _closuresId;
+ // face integration point in the coordinate of the left and right element (one fullMatrix per closure)
+ std::vector<fullMatrix<double> > _integrationPoints;
+ // XYZ gradient of the shape functions of both elements on the integrations points of the face
+ // (iQP*3+iXYZ , iFace*NPsi+iPsi)
+ fullMatrix<double> _dPsiDxOnQP;
+ const polynomialBasis *_fs;
+ const dgGroupOfElements &_elementGroup;
+ const dgGroupOfFaces &_faceGroup;
+ std::vector<int>_elementId;
+ // normals at integration points (N*Ni) x 3
+ fullMatrix<double> _normals;
+ public:
+ void addElement(int iElement, int iClosure);
+ dgGroupOfConnections(const dgGroupOfElements &elementGroup, const dgGroupOfFaces &face, int pOrder);
+ inline const polynomialBasis *getFunctionSpace() {return _fs;}
+ inline const dgGroupOfElements &getGroupOfElements() {return _elementGroup;}
+ inline int getElementId(int i) {return _elementId[i];}
+ inline MElement *getElement(int i) {return _elementGroup.getElement(_elementId[i]);}
+ inline const std::vector<int>& getClosure(int i) {return _closures[_closuresId[i]];}
+ inline fullMatrix<double>& getIntegrationPointsOnElement(int i) {return _integrationPoints[_closuresId[i]];}
+ inline fullMatrix<double>& getDPsiDxOnQP() {return _dPsiDxOnQP;}
+ inline fullMatrix<double>& getNormals() {return _normals;}
+ void init();
+};
+
class dgGroupOfFaces {
+ std::vector<dgGroupOfConnections*> _connections;
// normals always point outside left to right
// only used if this is a group of boundary faces
std::string _boundaryTag;
- const dgGroupOfElements &_groupLeft,&_groupRight;
- void addFace(const MFace &topoFace, int iElLeft, int iElRight);
- void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
- void addVertex(MVertex *topoVertex, int iElLeft, int iElRight);
// Two polynomialBases for left and right elements
// the group has always the same types for left and right
- const polynomialBasis *_fsLeft,*_fsRight, *_fsFace;
+ const polynomialBasis *_fsFace;
// N elements in the group
- std::vector<int>_left, _right;
std::vector<MElement *>_faces;
// Ni integration points, matrix of size Ni x 3 (u,v,weight)
fullMatrix<double> *_integration;
@@ -144,19 +170,6 @@ class dgGroupOfFaces {
// is characterized by a single integer which is the combination
// this closure is for the interpolation that MAY BE DIFFERENT THAN THE
// GEOMETRICAL CLOSURE !!!
- std::vector<std::vector<int> > _closuresLeft;
- std::vector<std::vector<int> > _closuresRight;
- std::vector<int> _closuresIdLeft;
- std::vector<int> _closuresIdRight;
- // face integration point in the coordinate of the left and right element (one fullMatrix per closure)
- std::vector<fullMatrix<double> > _integrationPointsLeft;
- std::vector<fullMatrix<double> > _integrationPointsRight;
- // XYZ gradient of the shape functions of both elements on the integrations points of the face
- // (iQP*3+iXYZ , iFace*NPsi+iPsi)
- fullMatrix<double> *_dPsiLeftDxOnQP;
- fullMatrix<double> *_dPsiRightDxOnQP;
- // normals at integration points (N*Ni) x 3
- fullMatrix<double> *_normals;
// detJac at integration points (N*Ni) x 1
fullMatrix<double> *_detJac;
// collocation matrices \psi_i (GP_j)
@@ -169,21 +182,6 @@ class dgGroupOfFaces {
//common part of the 3 constructors
void init(int pOrder);
public:
- inline const dgGroupOfElements &getGroupLeft()const {return _groupLeft; }
- inline const dgGroupOfElements &getGroupRight()const {return _groupRight; }
- inline int getElementLeftId (int i) const {return _left[i];};
- inline int getElementRightId (int i) const {return _right[i];};
- inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
- inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
- inline double getElementVolumeLeft(int iFace) const {return _groupLeft.getElementVolume(_left[iFace]);}
- inline double getElementVolumeRight(int iFace) const {return _groupRight.getElementVolume(_right[iFace]);}
- inline MElement* getFace (int iElement) const {return _faces[iElement];}
- inline const std::vector<int> &getClosureLeft(int iFace) const{ return _closuresLeft[_closuresIdLeft[iFace]]; }
- inline const std::vector<int> &getClosureRight(int iFace) const{ return _closuresRight[_closuresIdRight[iFace]];}
- inline fullMatrix<double> &getIntegrationOnElementLeft(int iFace) { return _integrationPointsLeft[_closuresIdLeft[iFace]];}
- inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _integrationPointsRight[_closuresIdRight[iFace]];}
-
- inline fullMatrix<double> &getNormals () const {return *_normals;}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder, int numVertices = -1);
dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder, int numVertices = -1);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
@@ -192,8 +190,6 @@ public:
virtual ~dgGroupOfFaces ();
inline bool isBoundary() const {return !_boundaryTag.empty();}
inline const std::string getBoundaryTag() const {return _boundaryTag;}
- inline fullMatrix<double> & getDPsiLeftDxMatrix() const { return *_dPsiLeftDxOnQP;}
- inline fullMatrix<double> & getDPsiRightDxMatrix() const { return *_dPsiRightDxOnQP;}
//this part is common with dgGroupOfElements, we should try polymorphism
inline int getNbElements() const {return _faces.size();}
inline int getNbNodes() const {return _collocation->size2();}
@@ -203,12 +199,34 @@ public:
inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iGaussPoint,iElement);}
inline double getInterfaceSurface (int iFace)const {return (*_interfaceSurface)(iFace,0);}
+ const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
+ inline MElement* getFace (int iElement) const {return _faces[iElement];}
+ // duplicate
+private:
+ void addFace(const MFace &topoFace, int iElLeft, int iElRight);
+ void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
+ void addVertex(MVertex *topoVertex, int iElLeft, int iElRight);
+public:
//keep this outside the Algorithm because this is the only place where data overlap
+ inline fullMatrix<double> &getNormals () const {return _connections[0]->getNormals();}
void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
void mapLeftFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft);
void mapRightFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vRight);
- const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
+ inline fullMatrix<double> & getDPsiLeftDxMatrix() const { return _connections[0]->getDPsiDxOnQP();}
+ inline fullMatrix<double> & getDPsiRightDxMatrix() const { return _connections[1]->getDPsiDxOnQP();}
+ inline const dgGroupOfElements &getGroupLeft()const {return _connections[0]->getGroupOfElements(); }
+ inline const dgGroupOfElements &getGroupRight()const {return _connections[1]->getGroupOfElements(); }
+ inline int getElementLeftId (int i) const {return _connections[0]->getElementId(i);}
+ inline int getElementRightId (int i) const {return _connections[1]->getElementId(i);}
+ inline MElement* getElementLeft (int i) const {return _connections[0]->getElement(i);}
+ inline MElement* getElementRight (int i) const {return _connections[1]->getElement(i);}
+ inline double getElementVolumeLeft(int iFace) const {return _connections[0]->getGroupOfElements().getElementVolume(getElementLeftId(iFace));}
+ inline double getElementVolumeRight(int iFace) const {return _connections[1]->getGroupOfElements().getElementVolume(getElementRightId(iFace));}
+ inline const std::vector<int> &getClosureLeft(int iFace) const{ return _connections[0]->getClosure(iFace);}
+ inline const std::vector<int> &getClosureRight(int iFace) const{ return _connections[1]->getClosure(iFace);}
+ inline fullMatrix<double> &getIntegrationOnElementLeft(int iFace) { return _connections[0]->getIntegrationPointsOnElement(iFace);}
+ inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _connections[1]->getIntegrationPointsOnElement(iFace);}
};
class dgGroupCollection {