dg: create new group collection by physical tag

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())

Solver/TESTCASES/Projection2D3D.lua

+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');
+

Solver/TESTCASES/prisms.geo

+//*********** 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]};
+

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);

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();
 };

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());