diff --git a/Solver/TESTCASES/implicit.lua b/Solver/TESTCASES/implicit.lua
new file mode 100644
index 0000000000000000000000000000000000000000..fa980b9d6afb0c9bf61064898ce304905175598c
--- /dev/null
+++ b/Solver/TESTCASES/implicit.lua
@@ -0,0 +1,16 @@
+model = GModel()
+model:load ('square2.msh')
+order = 1
+dimension = 2
+
+law = dgConservationLawAdvectionDiffusion:diffusionLaw(functionConstant({1}))
+law:addBoundaryCondition('Border',law:new0FluxBoundary())
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+solution = dgDofContainer(groups, law:getNbFields())
+residual = dgDofContainer(groups, law:getNbFields())
+r = dgResidualVolume (law)
+g = groups:getElementGroup(0)
+j = fullMatrix(3,3*g:getNbElements())
+r:compute1GroupWithJacobian(g,solution:getGroupProxy(0),residual:getGroupProxy(0),j)
+j:print("test")
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index b4144ca69fdb70e444c63b0f5b299d66b8da653b..79b9942a39059f704e5db5f7c2a8e21090e96bbe 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -121,7 +121,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
_PsiDPsiDXi = fullMatrix<double> (_dimUVW*nbQP, nbNodes*nbNodes);
//reditribution of the diffusive jacobian dimUVW*dimUVW*nbIntegrationPoints x nbNodes*nbNodes
- _dPsiDXDPsiDXi = fullMatrix<double> (_dimUVW*_dimUVW*nbQP, nbNodes*nbNodes);
+ _dPsiDXDPsiDXi = fullMatrix<double> (nbNodes*nbNodes, _dimUVW*_dimUVW*nbQP);
for (int xi=0;xi<nbQP; xi++) {
_fs.df((*_integration)(xi,0),
@@ -145,7 +145,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
}
for (int alpha=0; alpha<_dimUVW; alpha++) for (int beta=0; beta<_dimUVW; beta++) {
for (int i=0; i<nbNodes; i++) for (int j=0; j<nbNodes; j++) {
- _dPsiDXDPsiDXi((alpha*_dimUVW+beta)*nbQP+xi, i*nbNodes+j) = g[i][alpha]*g[j][beta]*weight;
+ _dPsiDXDPsiDXi(i*nbNodes+j, (alpha*_dimUVW+beta)*nbQP+xi) = g[i][alpha]*g[j][beta]*weight;
}
}
for (int alpha=0; alpha<_dimUVW; alpha++){
@@ -1212,6 +1212,8 @@ void dgGroupCollection::registerBindings(binding *b)
methodBinding *cm;
cb = b->addClass<dgGroupOfElements>("dgGroupOfElements");
cb->setDescription("a group of element sharing the same properties");
+ cm = cb->addMethod("getNbElements",&dgGroupOfElements::getNbElements);
+ cm->setDescription("return the number of elements in the group");
cb = b->addClass<dgGroupOfFaces>("dgGroupOfFaces");
cb->setDescription("a group of faces. All faces of this group are either interior of the same group of elements or at the boundary between two group of elements");
cb = b->addClass<dgGroupCollection>("dgGroupCollection");
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index e64dc9678de989f1e67f981c81a16046a6184cd4..3ef7ab80e864a59e6ebf228063a0c8d415e93040 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -135,10 +135,8 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
fullMatrix<double> fuvwe;
fullMatrix<double> source;
fullMatrix<double> dPsiDx, dofs, gradSolProxy;
- int nColA = group.getDimUVW()*group.getNbIntegrationPoints();
- int nColB = group.getDimUVW()*group.getDimUVW()*group.getNbIntegrationPoints();
- fullMatrix<double> A (_nbFields*_nbFields, nColA*group.getNbElements());
- fullMatrix<double> B (_nbFields*_nbFields, nColB*group.getNbElements());
+ fullMatrix<double> A (group.getDimUVW()*group.getNbIntegrationPoints(), _nbFields*_nbFields*group.getNbElements());
+ fullMatrix<double> B (group.getDimUVW()*group.getDimUVW()*group.getNbIntegrationPoints(), _nbFields*_nbFields*group.getNbElements());
fullMatrix<double> a, b;
functionDerivator *dDiffusiveFluxdGradU,*dConvectiveFluxdU;
if(_diffusiveFlux)
@@ -149,8 +147,8 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
_solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
- a.setAsProxy(A, iElement*nColA, nColA);
- b.setAsProxy(B, iElement*nColB, nColB);
+ a.setAsProxy(A, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
+ b.setAsProxy(B, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
if(_gradientSolutionQPe.somethingDependOnMe()){
dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
@@ -194,27 +192,26 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
}
}
int nQP = group.getNbIntegrationPoints();
- int dim = group.getDimUVW();
+ int dimUVW = group.getDimUVW();
+ int dimXYZ = group.getDimXYZ();
+ b.scale(0);
if(_diffusiveFlux) {
- printf("ok %i\n", __LINE__);
dDiffusiveFluxdGradU->compute();
- printf("ok %i\n", __LINE__);
- exit(0);
- for (int alpha=0; alpha < dim; alpha++) for (int beta=0; beta < dim; beta++) {
- for(int x=0;x <group.getDimXYZ();x++) for(int y=0;y <group.getDimXYZ();x++) {
+ for (int alpha=0; alpha < dimUVW; alpha++) for (int beta=0; beta < dimUVW; beta++) {
+ for(int x=0;x <dimXYZ;x++) for(int y=0;y <dimXYZ;y++) {
for (int xi=0; xi <group.getNbIntegrationPoints(); xi++) {
const double invJx = group.getInvJ (iElement, xi, alpha, x);
const double invJy = group.getInvJ (iElement, xi, beta, y);
const double detJ = group.getDetJ (iElement, xi);
const double factor = invJx * invJy * detJ;
for (int k=0; k< _nbFields; k++) for (int l=0; l<_nbFields; l++) {
- b(k*_nbFields+l,(alpha*dim+beta)*nQP+xi) = dDiffusiveFluxdGradU->get(k*dim+alpha,l*dim+beta,xi)*factor;
+ b((alpha*dimUVW+beta)*nQP+xi,k*_nbFields+l) += dDiffusiveFluxdGradU->get(xi, k*dimXYZ+x, l*dimXYZ+y)*factor;
}
}
}
}
}
- if(_convectiveFlux) {
+ /*if(_convectiveFlux) {
dConvectiveFluxdU->compute();
for (int alpha=0; alpha < dim; alpha++) {
for(int x=0;x <group.getDimXYZ();x++) {
@@ -228,33 +225,33 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
}
}
}
- }
+ }*/
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
- if(_convectiveFlux || _diffusiveFlux){
+ /*if(_convectiveFlux || _diffusiveFlux){
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
}
}
if(_sourceTerm){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
- }
+ }*/
int nbNodes = group.getNbNodes();
- fullMatrix<double> jacobianK (_nbFields*_nbFields,nbNodes*nbNodes);
- if (_convectiveFlux) {
- jacobianK.gemm(group.getPsiDPsiDXi(),B);
- }
- if (_convectiveFlux) {
- jacobianK.gemm(group.getDPsiDXDPsiDXi(),A);
+ fullMatrix<double> jacobianK (nbNodes*nbNodes,_nbFields*_nbFields*group.getNbElements());
+ /*if (_convectiveFlux) {
+ jacobianK.gemm(group.getPsiDPsiDXi(),A);
+ }*/
+ if (_diffusiveFlux) {
+ jacobianK.gemm(group.getDPsiDXDPsiDXi(),B);
}
fullMatrix<double> jacobianE, jacobianKE;
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
jacobianKE.setAsProxy(jacobianK, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
- jacobianKE.setAsProxy(jacobian, iElement*_nbFields*nbNodes, _nbFields*nbNodes);
+ jacobianE.setAsProxy(jacobian, iElement*_nbFields*nbNodes, _nbFields*nbNodes);
for (int k=0; k<_nbFields;k++) for (int l=0;l<_nbFields;l++) {
for(int i=0; i<nbNodes; i++) for(int j=0; j<nbNodes; j++) {
- jacobianE(l*nbNodes+j, k*nbNodes+i) = jacobianKE(k*_nbFields+l, i*nbNodes+j);
+ jacobianE(l*nbNodes+j, k*nbNodes+i) = jacobianKE(i*nbNodes+j, k*_nbFields+l);
}
}
}
diff --git a/Solver/functionDerivator.cpp b/Solver/functionDerivator.cpp
index b64897e65b2223ce61d0265823104f4a3b87816a..88d0e93541396f3f732d21b4558198285992172a 100644
--- a/Solver/functionDerivator.cpp
+++ b/Solver/functionDerivator.cpp
@@ -4,7 +4,6 @@
void functionDerivator::compute() {
_xRef = _x();
_fRef = _f();
- printf("f.size = %i %i\n_x.size = %i %i\n",_f().size1(), _f().size2(), _x().size1(), _x().size2());
_dfdx = fullMatrix<double>(_f().size1(),_f().size2()*_x().size2());
for (int j=0;j<_xRef.size2();j++) {
_xDx = _xRef;
@@ -16,5 +15,4 @@ void functionDerivator::compute() {
_dfdx(i, k*_xRef.size2()+j) = (_f(i,k)-_fRef(i,k))/_epsilon;
}
_x.set(_xRef);
- _dfdx.print();
};