Commit eec68ea7 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

introduce FEElement and FEFields

parent b40fadaa
......@@ -23,6 +23,7 @@
set(SRC
fluid_problem.c
fluid_problem_concentration.c
fluid_problem_fem.c
mesh.c
mesh_find.c
../tools/quadtree.c
......
......@@ -30,11 +30,9 @@
#include "mesh_find.h"
#include "fluid_problem_fem.h"
#define P D
#define P D
#define U 0
#define LOCAL_MATRIX(i,j,a,b) local_matrix[((i)*n_fields*N_SF+N_SF*(a)+(j))*n_fields+(b)]
static MeshTree *get_mesh_tree(FluidProblem *problem)
{
if (problem->mesh_tree == NULL)
......@@ -42,7 +40,7 @@ static MeshTree *get_mesh_tree(FluidProblem *problem)
return problem->mesh_tree;
}
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip);
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double dsol[][D], const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip);
void fluid_problem_interpolate_rho_and_nu(const FluidProblem *problem, double a, double *rho, double *mu) {
if (problem->n_fluids == 1){
*rho = problem->rho[0];
......@@ -83,18 +81,12 @@ inline static double node_volume_from_detj(double detj) {
int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}
void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *porosity = problem->porosity;
const double *solution = problem->solution;
int drag_in_stab = problem->drag_in_stab;
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
double *s = malloc(sizeof(double)*n_fields);
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
const FEFields *fields = problem->fields;
for (int ip=0; ip < problem->n_particles; ++ip) {
int iel = problem->particle_element_id[ip];
if(iel < 0){
......@@ -103,42 +95,30 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
continue;
}
double *xi = problem->particle_uvw + D*ip;
double phi[N_SF];
shape_functions(xi,phi);
const int *el = &mesh->elements[iel*N_N];
double dxidx[D][D], dphi[N_SF][D];
double sf[fields->local_size];
fe_fields_f(fields, xi, sf);
double sfporosity[problem->field_porosity->local_size];
fe_fields_f(problem->field_porosity, xi, sfporosity);
double dxidx[D][D];
double detj = mesh_dxidx(mesh, iel, dxidx);
grad_shape_functions(dxidx, dphi);
double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
sold[i] = 0;
for (int j = 0; j < D; ++j) {
ds[i*D+j] = 0;
}
}
double c = 0;
double dsf[fields->local_size][D];
fe_fields_df(fields, xi, dxidx, dsf);
double dsfporosity[problem->field_porosity->local_size][D];
fe_fields_df(problem->field_porosity, xi, dxidx, dsfporosity);
double *local_vector = all_local_vector ? &all_local_vector[fields->local_size*iel] : NULL;
double *local_matrix = all_local_matrix ? &all_local_matrix[fields->local_size*fields->local_size*iel] : NULL;
double c, cold, dc[D];
double s[n_fields], ds[n_fields][D];
fe_fields_eval_grad(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
double sold[n_fields];
fe_fields_eval(fields, mesh, iel, sf,solution_old, sold);
fe_fields_eval_grad(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
fe_fields_eval(problem->field_porosity, mesh, iel, sfporosity, problem->oldporosity, &cold);
double a = 0;
double cold = 0;
double dc[D] = {0};
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
cold += problem->oldporosity[el[i]]*phi[i];
if (problem->n_fluids == 2) {
a += problem->concentration[iel*N_N+i]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
double dofold = solution_old[el[i]*n_fields+j];
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
for (int k = 0; k < D; ++k)
ds[j*D+k] += dphi[i][k]*dof;
}
for (int k=0; k<D; ++k){
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
if (problem->n_fluids == 2) {
double sfconcentration[problem->field_concentration->local_size];
fe_fields_f(problem->field_concentration, xi, sfconcentration);
fe_fields_eval(problem->field_concentration, mesh, iel, problem->concentration, sfconcentration, &a);
}
double f[D],dfdu,dfddp;
particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
......@@ -147,40 +127,51 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
double rho,nu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
double pspg = drag_in_stab?problem->taup[iel]/rho:0;
for (int iD = 0; iD < D; ++iD) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
local_vector[iphi+N_SF*P] += pspg*dphi[iphi][iD]*f[iD];
for (int jD = 0; jD < D; ++jD) {
double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][jD]*f[iD];
double f0[n_fields], f1[n_fields][D], f00[n_fields*n_fields],f11[n_fields*n_fields][D][D],f10[n_fields*n_fields][D],f01[n_fields*n_fields][D];
for (int i = 0; i < n_fields; ++i) {
f0[i] = 0;
for (int id = 0; id < D; ++id) {
f1[i][id] = 0;
}
for (int j = 0; j < n_fields; ++j) {
int r = i*n_fields+j;
f00[r] = 0;
for (int id = 0; id < D; ++id) {
f01[r][id] = 0;
f10[r][id] = 0;
for (int jd = 0; jd < D; ++jd) {
f11[r][id][jd] = 0;
}
}
}
}
for (int d = 0; d < D; ++d) {
f0[U+d] += f[d];
f1[P][d] += pspg*f[d];
for (int e = 0; e < D; ++e) {
double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
f1[U+d][d] += supg*f[d];
}
}
fe_fields_add_to_local_vector(problem->fields, f0, f1, sf, dsf, 1, local_vector);
if (!local_matrix)
continue;
for (int iD = 0; iD < D; ++iD){
for (int iphi = 0; iphi < N_SF; ++iphi){
for (int jphi = 0; jphi < N_SF; ++jphi){
LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += phi[iphi]*phi[jphi]*dfdu;
LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][iD]*dfddp;
LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu;
LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][iD]*dfddp;
for (int jD = 0; jD < D; ++jD) {
double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += supg*dphi[iphi][jD]*phi[jphi]*dfdu;
LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][jD]*dphi[jphi][iD]*dfddp;
}
}
for (int d = 0; d < D; ++d) {
f00[(U+d)*n_fields+U+d] = dfdu;
f01[(U+d)*n_fields+P][d] = dfddp;
f10[P*n_fields+U+d][d] = pspg*dfdu;
f11[P*n_fields+P][d][d] = pspg*dfddp;
for (int e = 0; e < D; ++e) {
double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
f10[(U+d)*n_fields+U+d][e] += supg*dfdu;
f11[(U+d)*n_fields+P][d][e] += supg*dfddp;
}
}
fe_fields_add_to_local_matrix(problem->fields, f00, f01, f10, f11, sf, dsf, 1, local_matrix);
}
free(s);
free(ds);
free(sold);
}
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double *f01)
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double ds[][D], const double sold[D], const double dsold[][D], const double *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double f01[][D])
{
const int vid = wbnd->vid;
const int pid = wbnd->pid;
......@@ -197,7 +188,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
double dcn = 0;
for (int iD = 0; iD < D; ++iD) {
u[iD] = s[U+iD];
dp[iD] = ds[P*D+iD];
dp[iD] = ds[P][iD];
unold += sold[U+iD]*n[iD];
unmesh += mesh_velocity[iD]*n[iD];
uext[iD] = (vid<0?s[U+iD]:data[vid+iD]*c);
......@@ -224,7 +215,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
for (int iD = 0; iD < D; ++iD) {
for (int jD = 0; jD < D; ++jD)
c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
c_du_o_c[iD][jD] = ds[U+iD][jD] -u[iD]/c*dc[jD];
}
if(wbnd->compute_viscous_term == 1){
for (int id = 0; id < D; ++id) {
......@@ -234,8 +225,8 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
f01[(U+id)*n_fields+U+id][jd] -= mu*n[jd];
f01[(U+id)*n_fields+U+jd][id] -= mu*n[jd];
}
}
}
......@@ -337,13 +328,13 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
}
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double dsol[][DIMENSION], const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
const double *contact = &problem->particle_contact[ip*D];
double u[D], uold[D], dp[D];
for (int iD = 0; iD < D; ++iD) {
u[iD] = sol[U+iD];
uold[iD] = solold[U+iD];
dp[iD] = dsol[P*D+iD];
dp[iD] = dsol[P][iD];
}
double mu, rho;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
......@@ -381,12 +372,25 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
*dfddp = gammastar*dt/mass*vol;//-vol;
}
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double dsol[][D], const double *solold, const double dsolold[][D], const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double f1[][D], double f00[], double f10[][D], double f01[][D], double f11[][D][D]) {
double p = sol[P];
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
const int n_fields = fluid_problem_n_fields(problem);
double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
for (int i = 0; i < n_fields; ++i) {
for (int j = 0; j < n_fields; ++j) {
int r = i*n_fields+j;
f00[r] = 0;
for (int id = 0; id < D; ++id) {
f01[r][id] = 0;
f10[r][id] = 0;
for (int jd = 0; jd < D; ++jd) {
f11[r][id][jd] = 0;
}
}
}
}
if(problem->stab_param != 0)
taup = tauc = 0;
double divu = 0;
......@@ -396,10 +400,10 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
uold[iD] = solold[U+iD];
nold += uold[iD]*uold[iD];
umesh[iD] = mesh_velocity[iD];
dp[iD] = dsol[P*D+iD];
dp[iD] = dsol[P][iD];
for (int jD = 0; jD < D; ++jD) {
du[iD][jD] = dsol[(U+iD)*D+jD];
duold[iD][jD] = dsolold[(U+iD)*D+jD];
du[iD][jD] = dsol[U+iD][jD];
duold[iD][jD] = dsolold[U+iD][jD];
}
divu += du[iD][iD];
}
......@@ -416,12 +420,12 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
- bf[i]*c
+ drag*u[i];
//+ 5.3e5*u[i];
f00[(U+i)*n_fields+U+i] = drag;//5.3e5;
f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
f01[(U+i)*n_fields+P][i] += c;
if(problem->temporal) {
f0[U+i] += rho*(u[i]-uold[i])/dt;
f00[(U+i)*n_fields+U+i] += rho/dt;
}
f01[(U+i)*n_fields*D+P*D+i] = c;
if (problem->advection) {
// advection :
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
......@@ -430,58 +434,49 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
f0[U+i] -= rho/c*(umesh[j]*du[i][j]);
f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j];
f01[(U+i)*n_fields*D+(U+i)*D+j] -= rho/c*umesh[j];
f01[(U+i)*n_fields+(U+i)][j] += rho/c*uold[j];
f01[(U+i)*n_fields+(U+i)][j] -= rho/c*umesh[j];
}
}
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
f10[((U+i)*D+j)*n_fields+U+j] += -mu*dc[i]/c;
f10[((U+i)*D+j)*n_fields+U+i] += -mu*dc[j]/c;
f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu;
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+j] += mu;
f1[U+i][j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
f10[(U+i)*n_fields+U+j][j] += -mu*dc[i]/c;
f10[(U+i)*n_fields+U+i][j] += -mu*dc[j]/c;
f11[(U+i)*n_fields+U+j][j][i] += mu;
f11[(U+i)*n_fields+U+i][j][j] += mu;
// SUPG
double supg = uold[j]*taup;
f1[(U+i)*D+j] += supg*f0[U+i];
f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i];
f1[U+i][j] += supg*f0[U+i];
f10[(U+i)*n_fields+U+i][j] += supg*f00[(U+i)*n_fields+U+i];
f11[(U+i)*n_fields+P][j][i] += supg*f01[(U+i)*n_fields+P][i];
for (int k = 0; k < D; ++k)
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
f11[(U+i)*n_fields+(U+i)][j][k] += supg*f01[(U+i)*n_fields+(U+i)][k];
}
double lsic = tauc*rho;
f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
f1[U+i][i] += ((c-cold)/dt+divu)*lsic;
for (int j = 0; j < D; ++j) {
f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
f11[(U+i)*n_fields+(U+j)][i][j] += lsic;
}
f1[P*D+i] = -u[i];
f10[(P*D+i)*n_fields+U+i] += -1;
f1[P][i] = -u[i];
f10[P*n_fields+U+i][i] += -1;
// PSPG
double pspg = taup/rho;
f1[P*D+i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
f11[(P*D+i)*(n_fields*D)+P*D+i] += pspg*(f01[(U+i)*n_fields*D+P*D+i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
f10[(P*D+i)*n_fields+U+i] += pspg*f00[(U+i)*n_fields+(U+i)];
f1[P][i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
f11[P*n_fields+P][i][i] += pspg*(f01[(U+i)*n_fields+P][i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
f10[P*n_fields+U][i] += pspg*f00[(U+i)*n_fields+(U+i)];
for (int j = 0; j < D; ++j) {
f11[(P*D+i)*n_fields*D+(U+i)*D+j] = pspg*f01[(U+i)*n_fields*D+(U+i)*D+j];
f11[P*n_fields+U+i][i][j] = pspg*f01[(U+i)*n_fields+(U+i)][j];
}
}
}
void weak_boundary_values(const FEFields *fields, const Mesh *mesh, const MeshBoundary *bnd, const WeakBoundary *wbnd, double *data);
static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
const int n_fields = fluid_problem_n_fields(problem);
const size_t local_size = N_SF*n_fields;
double *s = malloc(sizeof(double)*n_fields);
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
double *dsold = malloc(sizeof(double)*n_fields*D);
double *f0 = malloc(sizeof(double)*n_fields);
double *f00 = malloc(sizeof(double)*n_fields*n_fields);
double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
double i_bnd = 0;
double iv_bnd_up = 0;
double iv_bnd_down = 0;
const size_t local_size = problem->fields->local_size;
const FEFields *fields = problem->fields;
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
int bndid= -1;
......@@ -496,130 +491,74 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
MeshBoundary *bnd = &problem->boundaries[bndid];
int n_value = weak_boundary_n_values(wbnd);
double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
weak_boundary_values(mesh, bnd, wbnd, data);
weak_boundary_values(problem->fields, mesh, bnd, wbnd, data);
for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
const int iel = interface[0];
const int icl = interface[1];
const int *cl = elbnd[icl];
int nodes[N_LSF];
double *x[N_LSF];
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
const int *el = &mesh->elements[iel*N_N];
for (int i = 0; i < D; ++i){
nodes[i] = el[cl[i]];
x[i] = &mesh->x[nodes[i]*3];
}
double dxidx[D][D], dphi[N_N][D];
mesh_dxidx(mesh,iel,dxidx);
double n[D],detbnd;
get_normal_and_det(x,n,&detbnd);
grad_shape_functions(dxidx, dphi);
double dc[D] = {0};
for (int j = 0; j < n_fields*D; ++j) {
ds[j] = 0;
}
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
double dofold = solution_old[el[i]*n_fields+j];
for (int k = 0; k < D; ++k){
ds[j*D+k] += dphi[i][k]*dof;
dsold[j*D+k] += dphi[i][k]*dofold;
}
}
for (int k=0; k<D; ++k){
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
double n[D];
double detbnd = get_normal(mesh, interface, n);
const FEElement *element = mesh->element;
element->df(QP[0], dxidx, dphi);
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
double dp[D] = {0};
double umesh[D] = {0};
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
double phi[DIMENSION];
l_shape_functions(LQP[iqp],phi);
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
sold[i] = 0;
double xiel[D];
param_boundary_to_element(interface[1],LQP[iqp], xiel);
double s[fields->n], sold[fields->n], ds[fields->n][D], dsold[fields->n][D];
double sf[fields->local_size], dsf[fields->local_size][D];
fe_fields_f(fields, xiel, sf);
fe_fields_df(fields, QP[iqp], dxidx, dsf);
fe_fields_eval_grad(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
fe_fields_eval_grad(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
double dc[D], c, sfporosity[problem->field_porosity->local_size], dsfporosity[problem->field_porosity->local_size][D] ;
fe_fields_f(problem->field_porosity, QP[iqp], sfporosity);
fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
fe_fields_eval_grad(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
double f0[fields->n], f00[fields->n*fields->n],f01[fields->n*fields->n][D];
for (int i = 0; i < fields->n; ++i) {
f0[i] = 0;
for (int j = 0; j < fields->n; ++j) {
int r = i*fields->n+j;
f00[r] = 0;
for (int id = 0; id < D; ++id) {
f01[r][id] = 0;
}
}
}
for (int i = 0; i < n_fields*n_fields; ++i) {
f00[i] = 0;
}
for (int i = 0; i < n_fields*n_fields*D; ++i) {
f01[i] = 0;
}
double c = 0;
double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
double a = 0;
for (int i = 0; i < DIMENSION; ++i) {
c += problem->porosity[nodes[i]]*phi[i];
if (problem->n_fluids == 2) {
a += problem->concentration[iel*N_N+cl[i]]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[nodes[i]*n_fields+j];
double dofold = solution_old[nodes[i]*n_fields+j];
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
}
for(int j = 0; j < D; ++j){
double dofmesh = problem->mesh_velocity[nodes[i]*D+j];
umesh[j] += phi[i]*dofmesh;
if (problem->n_fluids==2) {
double sfconcentration[problem->field_concentration->local_size];
fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
fe_fields_eval(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
}
const FEElement *mesh_element = problem->mesh->element;
double meshf[mesh_element->nlocal];
mesh_element->f(QP[iqp], meshf);
double umesh[D] = {0};
for (int i = 0; i < mesh_element->nlocal; ++i) {
for (int j=0; j<D; ++j) {
umesh[j] += meshf[i]*problem->mesh_velocity[el[i]*D+j];
}
}
double rho, mu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
const double jw = LQW[iqp]*detbnd;
if (wbnd->type != BND_WALL){
for (int i = 0; i < D; ++i){
if (wbnd->vid<0) {
i_bnd -= a*s[U+i]*n[i]*jw*dt;
iv_bnd_up -= s[U+i]*n[i]*jw;
}
else {
iv_bnd_down -= dataqp[wbnd->vid+i]*n[i]*jw;
i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
}
}
}
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,umesh,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < D; ++iphi){
if (ifield<D){
problem->boundary_force[bndid*D + ifield] += phi[iphi]*f0[ifield]*jw;
}
local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
}
}
for (int jfield = 0; jfield < n_fields; ++jfield) {
for (int ifield = 0; ifield < n_fields; ++ifield){
for (int iphi = 0; iphi < N_LSF; ++iphi){
for (int jphi = 0; jphi < N_LSF; ++jphi){
double d = phi[jphi]*f00[ifield*n_fields+jfield];
LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*d;
}
for (int jphi = 0; jphi < N_SF; ++jphi){
double d = 0;
for (int jD = 0; jD < D; ++jD) {
d += dphi[jphi][jD]*f01[ifield*n_fields*D+jfield*D+jD];
}
LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += jw*phi[iphi]*d;
}
}
}
for (int ifield = 0; ifield < D; ++ifield) {
problem->boundary_force[bndid*D + ifield] += f0[ifield]*jw;
}
fe_fields_add_to_local_vector(fields, f0, NULL, sf, NULL, jw, local_vector);
fe_fields_add_to_local_matrix(fields, f00, f01, NULL, NULL, sf, dsf, jw, local_matrix);
}
}
free(data);
}
free(s);
free(ds);
free(sold);
free(dsold);
free(f0);