double v1dotn = 0.;
double u2dotn = 0.;
double v2dotn = 0.;
-
+
+ double ngradu1n = 0.;
+ double ngradv1n = 0.;
+ double ngradu2n = 0.;
+ double ngradv2n = 0.;
+
for (unsigned int d=0; d<dim; ++d)
{
u1dotn += n(d)*fe1.shape_value_component(j,k,d);
v1dotn += n(d)*fe1.shape_value_component(i,k,d);
u2dotn += n(d)*fe2.shape_value_component(j,k,d);
v2dotn += n(d)*fe2.shape_value_component(i,k,d);
+
+ ngradu1n += n*fe1.shape_grad_component(j,k,d)*n(d);
+ ngradv1n += n*fe1.shape_grad_component(i,k,d)*n(d);
+ ngradu2n += n*fe2.shape_grad_component(j,k,d)*n(d);
+ ngradv2n += n*fe2.shape_grad_component(i,k,d)*n(d);
}
for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
{
+ const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n(d);
+ const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n(d);
- // multiply by
- const double vi = fe1.shape_value_component(i,k,d)*(1-v1dotn);
- const double dnvi = n * fe1.shape_grad_component(i,k,d)*(1-v1dotn);
- const double ve = fe2.shape_value_component(i,k,d)*(1-v2dotn);
- const double dnve = n * fe2.shape_grad_component(i,k,d)*(1-v2dotn);
- const double ui = fe1.shape_value_component(j,k,d)*(1-u1dotn);
- const double dnui = n * fe1.shape_grad_component(j,k,d)*(1-u1dotn);
- const double ue = fe2.shape_value_component(j,k,d)*(1-u2dotn);
- const double dnue = n * fe2.shape_grad_component(j,k,d)*(1-u2dotn);
- M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
+ const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n(d);
+ const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n(d);
+
+ const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n(d);
+ const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n(d);
+
+ const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n(d);
+ const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n(d);
+
+ M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);