Commit 233d4c8f authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

transpose local matrix to have block structure

parent a9146bb2
......@@ -121,10 +121,13 @@ HXTStatus hxtLinearSystemPETScAddToMatrix(HXTLinearSystemPETSc *system, int el0,
int *map1 = malloc(sizeof(int)*nn*nf);
for(int i=0; i<nn; ++i){
for (int j=0; j<nf; ++j){
map0[j*nn+i] = system->nodeMap[e0[i]]*nf+j;
map1[j*nn+i] = system->nodeMap[e1[i]]*nf+j;
//map0[j*nn+i] = system->nodeMap[e0[i]]*nf+j;
//map1[j*nn+i] = system->nodeMap[e1[i]]*nf+j;
map0[i*nf+j] = system->nodeMap[e0[i]]*nf+j;
map1[i*nf+j] = system->nodeMap[e1[i]]*nf+j;
}
}
//HXT_PETSC_CHECK(MatSetValues(system->a,nn*nf,map0,nn*nf,map0,localMatrix,ADD_VALUES));
HXT_PETSC_CHECK(MatSetValues(system->a,nn*nf,map0,nn*nf,map1,localMatrix,ADD_VALUES));
free(map0);
free(map1);
......
......@@ -410,9 +410,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int jphi = 0; jphi < N_SF; ++jphi) {
double phij = phi[jphi];
const double *dphij = dphi[jphi];
int P = N_SF*DIMENSION;
int N2 = N_SF*N_SF;
int IDX = iphi*N2+jphi;
int P = DIMENSION;
//int N2 = N_SF*N_SF;
//int IDX = iphi*N2+jphi;
int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double dphiidphij = 0;
double dtau[DIMENSION];
......@@ -425,11 +427,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
dphiidtau += dphii[j]*dtau[j];
}
for (int j = 0; j < DIMENSION; ++j) {
int U = N_SF*j;
int U = j;
LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
LOCAL_MATRIX(P,U) += jw*phii*dphij[j];
for (int k = 0; k < DIMENSION; ++k) {
int V = N_SF*k;
int V = k;
LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
}
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
......
......@@ -76,7 +76,7 @@ fluid.export(outputdir,0,0)
ii = 0
tic = time.clock()
#forces = g*p.mass()
while t < tEnd :
while t < tEnd and ii < 50:
ticf = time.clock()
forces = fluid.compute_node_force(dt)
tocf = time.clock()
......
......@@ -76,7 +76,7 @@ print("RHOP = %g" % rhop)
outf = 50
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Top",0,0.),("Box",1,0.),("Box",0,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,autoEpsilon,strong_boundaries)
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export(outputdir,0,0)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment