From c86857e6f1b6768c5ecf40e1f7aade70f8542f3d Mon Sep 17 00:00:00 2001
From: Tuomas Karna <tuomas.karna@uclouvain.be>
Date: Wed, 10 Mar 2010 09:13:25 +0000
Subject: [PATCH] dg: create new group collection by physical tag
---
Solver/TESTCASES/Diffusion3D.lua | 16 +++----------
Solver/TESTCASES/Projection2D3D.lua | 34 +++++++++++++++++++++++++++
Solver/TESTCASES/prisms.geo | 36 +++++++++++++++++++++++++++++
Solver/dgGroupOfElements.cpp | 32 ++++++++++++++++++++++---
Solver/dgGroupOfElements.h | 4 +++-
Solver/dgSystemOfEquations.cpp | 2 +-
6 files changed, 106 insertions(+), 18 deletions(-)
create mode 100644 Solver/TESTCASES/Projection2D3D.lua
create mode 100644 Solver/TESTCASES/prisms.geo
diff --git a/Solver/TESTCASES/Diffusion3D.lua b/Solver/TESTCASES/Diffusion3D.lua
index ac7550ef65..fb8f894bf0 100644
--- a/Solver/TESTCASES/Diffusion3D.lua
+++ b/Solver/TESTCASES/Diffusion3D.lua
@@ -21,23 +21,13 @@ end
-- conservation law
-- advection speed
-v=fullMatrix(3,1);
-v:set(0,0,0)
-v:set(1,0,0)
-v:set(2,0,0)
-
-nu=fullMatrix(1,1);
-nu:set(0,0,0.01)
-
+vFunc = functionConstant({0.1,0.1,0});
-- diffusion function
+nuFunc = functionConstant({0});
-law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),functionConstant(nu):getName())
+law = dgConservationLawAdvectionDiffusion(vFunc:getName(),nuFunc:getName())
dg:setConservationLaw(law)
--- boundary condition
-outside=fullMatrix(1,1)
-outside:set(0,0,0.15)
-
law:addBoundaryCondition('side1',law:new0FluxBoundary())
law:addBoundaryCondition('side2',law:new0FluxBoundary())
law:addBoundaryCondition('side3',law:new0FluxBoundary())
diff --git a/Solver/TESTCASES/Projection2D3D.lua b/Solver/TESTCASES/Projection2D3D.lua
new file mode 100644
index 0000000000..a62c38aba2
--- /dev/null
+++ b/Solver/TESTCASES/Projection2D3D.lua
@@ -0,0 +1,34 @@
+model = GModel()
+-- model:load('tetrahedra.geo')
+-- model:load('tetrahedra.msh')
+-- model:load('prisms.geo')
+-- model:load('prisms.msh')
+model:load('mixed.geo')
+model:load('mixed.msh')
+order=1
+
+-- initial condition
+function initial_condition( xyz , f )
+ for i=0,xyz:size1()-1 do
+ x = xyz:get(i,0)
+ y = xyz:get(i,1)
+ z = xyz:get(i,2)
+ f:set (i, 0, math.exp(-50*((x-0.4)^2+(y-0.3)^2)))
+-- f:set (i, 0, math.exp(-50*((x-0.4)^2+(z+0.3)^2)))
+ end
+end
+
+init_cond = functionLua(1, 'initial_condition', {'XYZ'}):getName()
+gc=dgGroupCollection(model,3,order)
+solTmp=dgDofContainer(gc,1)
+solTmp:L2Projection(init_cond)
+gc:buildGroupsOfInterfaces(model,2,order)
+solution=dgDofContainer(gc,1)
+solution:L2Projection(init_cond)
+solution:exportMsh('sol3d.msh');
+
+gc2D = gc:newByTag(model,2,order,{"top","side1"});
+solution2d=dgDofContainer(gc2D,1);
+solution2d:L2Projection(init_cond);
+solution2d:exportMsh('sol2d.msh');
+
diff --git a/Solver/TESTCASES/prisms.geo b/Solver/TESTCASES/prisms.geo
new file mode 100644
index 0000000000..9a191a1edc
--- /dev/null
+++ b/Solver/TESTCASES/prisms.geo
@@ -0,0 +1,36 @@
+//*********** prisms.geo *************//
+C = 1;
+Lup = 1;
+L = 1.;
+lc = .1;
+
+Point(1) = {0.0, 0.0, -Lup, lc};
+Point(2) = {C , 0.0, -Lup, lc};
+Point(3) = {C , C , -Lup, lc};
+Point(4) = {0.0, C , -Lup, lc};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+
+outpri[]= Extrude {0,0,1.2*L}{ Surface{6}; Layers{Ceil(1.2*L/lc)};Recombine;};
+// outv[]= Extrude {0,0,0.5*L}{ Surface{7}; Layers{5};Recombine;};
+// Printf("top surface = %g", outpri[0]);
+// Printf("volume = %g", outpri[1]);
+// Printf("side surfaces = %g %g %g %g", outpri[2], outpri[3], outpri[4], outpri[5]);
+
+Mesh.Algorithm3D=4; // frontal [lobes]
+
+Physical Surface("top") = {outpri[0]};
+Physical Surface("bottom") = {6};
+Physical Surface("side1") = {outpri[2]};
+Physical Surface("side2") = {outpri[3]};
+Physical Surface("side3") = {outpri[4]};
+Physical Surface("side4") = {outpri[5]};
+Physical Volume("volume") = {outpri[1]};
+
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index dfdb3fe8f0..7248b24dd5 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -332,7 +332,7 @@ void dgGroupOfConnections::init() {
// -) Elements of different types are separated
// -) Then those groups are split into partitions
-void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
+void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order, std::vector<std::string>* physicalTags)
{
if(_groupsOfElementsBuilt)
return;
@@ -345,7 +345,16 @@ void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
std::multimap<MElement*, short> &ghostsCells = _model->getGhostCells();
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *entity = entities[i];
- if(entity->dim() == dim){
+ bool correctPhysicalTag = false;
+ if(physicalTags==NULL || physicalTags->size()==0)
+ correctPhysicalTag = true;
+ else
+ for(unsigned int iPhy = 0; iPhy < entity->physicals.size(); iPhy++)
+ for(unsigned int jPhy = 0; jPhy < physicalTags->size(); jPhy++)
+ if (_model->getPhysicalName(dim,entity->physicals[iPhy]) == physicalTags->at(jPhy) )
+ correctPhysicalTag = true;
+
+ if(entity->dim() == dim && correctPhysicalTag){
for (int iel=0; iel<entity->getNumMeshElements(); iel++){
MElement *el=entity->getMeshElement(iel);
if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
@@ -1162,7 +1171,21 @@ dgGroupCollection::dgGroupCollection()
dgGroupCollection::dgGroupCollection(GModel *model, int dimension, int order)
{
_groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
- buildGroupsOfElements(model,dimension,order);
+ buildGroupsOfElements(model,dimension,order, NULL);
+}
+
+dgGroupCollection::dgGroupCollection(GModel *model, int dimension, int order, std::vector<std::string>* physicalTags)
+{
+ _groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
+ buildGroupsOfElements(model,dimension,order,physicalTags);
+}
+
+dgGroupCollection* dgGroupCollection::newByTag(GModel *model, int dimension, int order, std::vector<std::string> tags){
+ Msg::Info("Creating %dD groupCollection by tags.",dimension);
+ dgGroupCollection* subCollection = new dgGroupCollection(model,dimension,order,&tags);
+ if (subCollection->getNbElementGroups()==0)
+ Msg::Warning("groupCollection is empty.");
+ return subCollection;
}
dgGroupCollection::~dgGroupCollection()
@@ -1207,6 +1230,9 @@ void dgGroupCollection::registerBindings(binding *b)
cm = cb->addMethod("splitGroupsByVerticalLayer",&dgGroupCollection::splitGroupsByVerticalLayer);
cm->setDescription("Split the groups according vertical layer structure. The first is defined by the topLevelTags.");
cm->setArgNames("topLevelTags",NULL);
+ cm = cb->addMethod("newByTag",&dgGroupCollection::newByTag);
+ cm->setDescription("Creates a group collection of elements associated with given physical tag(s).");
+ cm->setArgNames("model","dimension","order","tags",NULL);
cm = cb->addMethod("getNbElementGroups", &dgGroupCollection::getNbElementGroups);
cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 0190ef5ccd..ed2ea0dc4d 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -244,7 +244,7 @@ class dgGroupCollection {
inline int getImageElementGroup(int partId, int i) const {return _elementsToSend[partId][i].first;}
inline int getImageElementPositionInGroup(int partId, int i) const {return _elementsToSend[partId][i].second;}
- void buildGroupsOfElements (GModel *model,int dimension, int order);
+ void buildGroupsOfElements (GModel *model,int dimension, int order, std::vector<std::string>* physicalTags);
void buildGroupsOfInterfaces ();
double splitGroupsForMultirate(int maxLevels,int bufferSize,dgConservationLaw *claw, dgDofContainer *solution);
@@ -253,7 +253,9 @@ class dgGroupCollection {
void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
dgGroupCollection(GModel *model, int dimension, int order);
+ dgGroupCollection(GModel *model, int dimension, int order, std::vector<std::string>* physicalTags);
dgGroupCollection();
+ static dgGroupCollection* newByTag(GModel* model, int dimension, int order, std::vector<std::string> tags);
static void registerBindings(binding *b);
~dgGroupCollection();
};
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 7f7d62219c..d0e4fa66bf 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -99,7 +99,7 @@ void dgSystemOfEquations::L2Projection (std::string functionName){
// ok, we can setup the groups and create solution vectors
void dgSystemOfEquations::setup(){
if (!_claw) throw;
- _groups.buildGroupsOfElements(_gm,_dimension,_order);
+ _groups.buildGroupsOfElements(_gm,_dimension,_order, NULL);
_groups.buildGroupsOfInterfaces();
_solution = new dgDofContainer(&_groups,_claw->getNbFields());
_rightHandSide = new dgDofContainer(&_groups,_claw->getNbFields());
--
GitLab