From cef8f7457a5c9602f4b704d3821bfa765aff6f2b Mon Sep 17 00:00:00 2001
From: Samuel Melchior <samuel.melchior@uclouvain.be>
Date: Fri, 12 Mar 2010 12:15:28 +0000
Subject: [PATCH] Inverse mass matrix directly multiply the sparse projection
matrix
---
Solver/TESTCASES/Projection2D3D.lua | 25 ++++++++++++--------
Solver/TESTCASES/prisms.geo | 6 ++---
Solver/dgDofContainer.cpp | 36 ++++++++++++++++++++++++++---
Solver/dgDofContainer.h | 1 +
4 files changed, 52 insertions(+), 16 deletions(-)
diff --git a/Solver/TESTCASES/Projection2D3D.lua b/Solver/TESTCASES/Projection2D3D.lua
index a62c38aba2..19b3954c30 100644
--- a/Solver/TESTCASES/Projection2D3D.lua
+++ b/Solver/TESTCASES/Projection2D3D.lua
@@ -1,10 +1,10 @@
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')
+model:load('prisms.geo')
+model:load('prisms.msh')
+-- model:load('mixed.geo')
+-- model:load('mixed.msh')
order=1
-- initial condition
@@ -13,11 +13,13 @@ function initial_condition( xyz , f )
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.5)^2+(y-0.5)^2)))
-- f:set (i, 0, math.exp(-50*((x-0.4)^2+(z+0.3)^2)))
end
end
+proj = linearSystemCSRGmmdouble()
+
init_cond = functionLua(1, 'initial_condition', {'XYZ'}):getName()
gc=dgGroupCollection(model,3,order)
solTmp=dgDofContainer(gc,1)
@@ -25,10 +27,13 @@ solTmp:L2Projection(init_cond)
gc:buildGroupsOfInterfaces(model,2,order)
solution=dgDofContainer(gc,1)
solution:L2Projection(init_cond)
-solution:exportMsh('sol3d.msh');
+solution:exportMsh('sol3d.msh')
+
+gc2D = gc:newByTag(model,2,order,{"top"})
+solution2d=dgDofContainer(gc2D,1)
+--solution2d:L2Projection(init_cond)
-gc2D = gc:newByTag(model,2,order,{"top","side1"});
-solution2d=dgDofContainer(gc2D,1);
-solution2d:L2Projection(init_cond);
-solution2d:exportMsh('sol2d.msh');
+solution2d:Mesh2Mesh_BuildL2Projection(proj,solution)
+solution2d:Mesh2Mesh_ApplyL2Projection(proj,solution)
+solution2d:exportMsh('sol2d.msh')
diff --git a/Solver/TESTCASES/prisms.geo b/Solver/TESTCASES/prisms.geo
index 9a191a1edc..2291c6120d 100644
--- a/Solver/TESTCASES/prisms.geo
+++ b/Solver/TESTCASES/prisms.geo
@@ -1,8 +1,8 @@
//*********** prisms.geo *************//
C = 1;
-Lup = 1;
L = 1.;
-lc = .1;
+Lup = L;
+lc = 0.1;
Point(1) = {0.0, 0.0, -Lup, lc};
Point(2) = {C , 0.0, -Lup, lc};
@@ -18,7 +18,7 @@ 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;};
+outpri[]= Extrude {0,0,L}{ Surface{6}; Layers{Ceil(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]);
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 78c2b4f387..6270552d01 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -304,9 +304,39 @@ void dgDofContainer::Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &pro
jGlobal += dGroup.getNbNodes();
}
}
+ multAddInverseMassMatrixL2Projection(projector);
// buffer.print();
}
+void dgDofContainer::multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<double> &projector){
+ dgGroupCollection* rGroups = this->getGroups();
+ int *startIndex;
+ int *columns;
+ double *values;
+ projector.getMatrix(startIndex,columns,values);
+
+ int iGlobal = 0;
+ for (int iGroup=0;iGroup<rGroups->getNbElementGroups();iGroup++) {// for 2d.groups
+ const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);
+ for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {// for elements
+ fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),startIndex[iGlobal+1]-startIndex[iGlobal]);
+ fullMatrix<double> buffer2 = fullMatrix<double> (rGroup.getNbNodes(),startIndex[iGlobal+1]-startIndex[iGlobal]);
+ fullMatrix<double> iMassEl;
+ for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
+ for (int i = startIndex[iGlobal+iNode]; i < startIndex[iGlobal+iNode+1] ; i++)
+ buffer(iNode,i-startIndex[iGlobal+iNode]) = values[i];
+ }
+ iMassEl.setAsProxy(rGroup.getInverseMassMatrix(),iElement*rGroup.getNbNodes(),rGroup.getNbNodes());
+ buffer2.gemm(iMassEl,buffer);
+ for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
+ for (int i = startIndex[iGlobal+iNode]; i < startIndex[iGlobal+iNode+1] ; i++)
+ values[i] = buffer2(iNode,i-startIndex[iGlobal+iNode]);
+ }
+ iGlobal+=rGroup.getNbNodes();
+ }
+ }
+}
+
void dgDofContainer::Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor){
scale(0.);
@@ -327,7 +357,7 @@ int nbFields = 1; // TO DEFINE !!!
for (int iGroup=0;iGroup<rGroups->getNbElementGroups();iGroup++) {// for 2d.groups
const dgGroupOfElements &rGroup = *rGroups->getElementGroup(iGroup);
- fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),rGroup.getNbElements() * _nbFields);
+// fullMatrix<double> buffer = fullMatrix<double> (rGroup.getNbNodes(),rGroup.getNbElements() * _nbFields);
this->getGroupProxy(iGroup).scale(0);
for (int iElement=0 ; iElement<rGroup.getNbElements() ;++iElement) {// for elements
for (int iNode = 0 ; iNode<rGroup.getNbNodes() ;++iNode) {
@@ -340,14 +370,14 @@ int nbFields = 1; // TO DEFINE !!!
int jElement = (jGlobal-dGroupsStartIGlobal[jGroup])/dGroups->getElementGroup(jGroup)->getNbNodes();
int jNode = jGlobal-dGroupsStartIGlobal[jGroup]-jElement*dGroups->getElementGroup(jGroup)->getNbNodes();
for (int m = 0; m < nbFields; m++){
- buffer(iNode,iElement*nbFields+m) += values[i]*(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m);
+ this->getGroupProxy(iGroup)(iNode,iElement*nbFields+m) += values[i]*(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m);
// printf("%d %d %d %f %f %f\n",iGlobal,jGlobal,i, buffer(iNode,iElement*nbFields+m),values[i],(donor.getGroupProxy(jGroup))(jNode,jElement*nbFields+m));
}
}
}
}//*/
// buffer.print();
- dgAlgorithm::multAddInverseMassMatrix(rGroup, buffer, this->getGroupProxy(iGroup));
+// dgAlgorithm::multAddInverseMassMatrix(rGroup, buffer, this->getGroupProxy(iGroup));
// this->getGroupProxy(iGroup).print();
}
}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 2565dd5cbf..e4531cef00 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -44,6 +44,7 @@ public:
void setAll(double v);
void L2Projection(std::string functionName);
void Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
+ void multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<double> &projector); // this method should be private
void Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
void exportMsh(const std::string name);
void exportGroupIdMsh();
--
GitLab