renamed some routines

src/feenox.h

@@ -2410,7 +2410,9 @@ extern double *feenox_fem_compute_x_at_gauss_and_update_var(element_t *e, unsign
 extern double feenox_fem_determinant(gsl_matrix *);
 extern gsl_matrix *feenox_fem_matrix_invert(gsl_matrix *direct, gsl_matrix *inverse);
 
-extern double feenox_fem_compute_w_det_at_gauss(element_t *e, unsigned int q, int integration);
+extern double feenox_fem_compute_w_det_at_gauss_integration(element_t *e, unsigned int q, int integration);
+#define feenox_fem_compute_w_det_at_gauss(e, q) feenox_fem_compute_w_det_at_gauss_integration((e), (q), feenox.pde.mesh->integration)
+
 extern gsl_matrix *feenox_fem_compute_C(element_t *e);
 extern gsl_matrix *feenox_fem_compute_H_c_at_gauss(element_t *e, unsigned int q, int integration);
 extern gsl_matrix *feenox_fem_compute_H_Gc_at_gauss(element_t *e, unsigned int q, int integration);
@@ -2425,7 +2427,9 @@ extern gsl_matrix *feenox_fem_compute_J_at_gauss_general(element_t *e, unsigned
 extern gsl_matrix *feenox_fem_compute_B_c(element_t *e, double *xi);
 extern gsl_matrix *feenox_fem_compute_B(element_t *e, double *xi);
 
-extern gsl_matrix *feenox_fem_compute_B_at_gauss(element_t *e, unsigned int q, int integration);
+extern gsl_matrix *feenox_fem_compute_B_at_gauss_integration(element_t *e, unsigned int q, int integration);
+#define feenox_fem_compute_B_at_gauss(e, q) feenox_fem_compute_B_at_gauss_integration((e), (q), feenox.pde.mesh->integration)
+
 extern gsl_matrix *feenox_fem_compute_B_Gc_at_gauss(element_t *e, unsigned int q, int integration);
 extern gsl_matrix *feenox_fem_compute_B_G_at_gauss(element_t *e, unsigned int q, int integration);
 
@@ -2478,7 +2482,8 @@ extern int feenox_problem_build_element_volumetric_gauss_point(element_t *, unsi
 extern int feenox_problem_build_elemental_objects_allocate(element_t *);
 extern int feenox_problem_build_elemental_objects_free(void);
 extern int feenox_problem_build_assemble(void);
-extern int feenox_problem_rhs_set(element_t *e, unsigned int q, double *value);
+extern int feenox_problem_rhs_add(element_t *e, unsigned int q, double *value);
+extern int feenox_problem_stiffness_add(element_t *e, unsigned int q, gsl_matrix *Kiq);
 
 // gradient.c
 extern int feenox_problem_gradient_compute(void);

src/mesh/integrate.c

@@ -103,7 +103,7 @@ double feenox_mesh_integral_over_element(element_t *this, mesh_t *mesh, function
   double integral = 0;
 
   for (unsigned int q = 0; q < this->type->gauss[mesh->integration].Q; q++) {
-    double wdet = feenox_fem_compute_w_det_at_gauss(this, q, mesh->integration);
+    double wdet = feenox_fem_compute_w_det_at_gauss_integration(this, q, mesh->integration);
     double val = 0;
     for (unsigned int j = 0; j < this->type->nodes; j++) {
       val += gsl_matrix_get(this->type->gauss[mesh->integration].H_c[q], 0, j) * feenox_vector_get(function->vector_value, this->node[j]->index_mesh);
@@ -157,7 +157,7 @@ double feenox_mesh_integral_function_node_gauss(function_t *function, mesh_t *me
     element_t *e = &mesh->element[i];
     if (is_element_local(mesh, e) && is_physical_group_right(e, physical_group)) {
       for (unsigned int q = 0; q < e->type->gauss[mesh->integration].Q; q++) {
-        double wdet = feenox_fem_compute_w_det_at_gauss(e, q, mesh->integration);
+        double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, mesh->integration);
 
         double xi = 0;
         for (unsigned int j = 0; j < e->type->nodes; j++) {
@@ -184,7 +184,7 @@ double feenox_mesh_integral_function_general_gauss(function_t *function, mesh_t
       for (unsigned int q = 0; q < e->type->gauss[mesh->integration].Q; q++) {
         // TODO: check if the integrand depends on space
         double *x = feenox_fem_compute_x_at_gauss(e, q, mesh->integration);
-        double wdet = feenox_fem_compute_w_det_at_gauss(e, q, mesh->integration);
+        double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, mesh->integration);
         integral += wdet * feenox_function_eval(function, x);
       }  
     }
@@ -246,7 +246,7 @@ double feenox_mesh_integral_expression_gauss(expr_t *expr, mesh_t *mesh, physica
           feenox_call(mesh_compute_normal(element));
         }  
 */
-        integral += feenox_fem_compute_w_det_at_gauss(e, q, mesh->integration) * feenox_expression_eval(expr);
+        integral += feenox_fem_compute_w_det_at_gauss_integration(e, q, mesh->integration) * feenox_expression_eval(expr);
       }
     }
   }

src/mesh/mesh.c

@@ -229,7 +229,7 @@ int feenox_instruction_mesh_read(void *arg) {
         for (size_t i = 0; i < physical_group->n_elements; i++) {
           element_t *element = &this->element[physical_group->element[i]];
           for (unsigned int q = 0; q < element->type->gauss[this->integration].Q; q++) {
-            double wdet = feenox_fem_compute_w_det_at_gauss(element, q, this->integration);
+            double wdet = feenox_fem_compute_w_det_at_gauss_integration(element, q, this->integration);
             gsl_matrix *H_c = feenox_fem_compute_H_c_at_gauss(element, q, this->integration);
 
             for (unsigned int j = 0; j < element->type->nodes; j++) {

src/mesh/physical_group.c

@@ -172,7 +172,7 @@ int feenox_physical_group_compute_volume(physical_group_t *this, const mesh_t *m
   for (size_t i = 0; i < this->n_elements; i++) {
     element_t *element = &mesh->element[this->element[i]];
     for (unsigned int q = 0; q < element->type->gauss[mesh->integration].Q; q++) {
-      double wdet = feenox_fem_compute_w_det_at_gauss(element, q, mesh->integration);
+      double wdet = feenox_fem_compute_w_det_at_gauss_integration(element, q, mesh->integration);
 
       for (size_t j = 0; j < element->type->nodes; j++) {
         double wh = wdet * gsl_matrix_get(element->type->gauss[mesh->integration].H_c[q], 0, j);

src/pdes/blas.c

@@ -55,7 +55,6 @@ int feenox_blas_BtB_accum(gsl_matrix *B, double alpha, gsl_matrix *R) {
 
 // R = alpha * B'*B
 int feenox_blas_BtB(gsl_matrix *B, double alpha, gsl_matrix *R) {
-//  printf("dgemm BtB\n");
   feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, alpha, B, B, 0.0, R));
   return FEENOX_OK;
 }

src/pdes/build.c

@@ -405,18 +405,29 @@ int feenox_problem_build_assemble(void) {
   return FEENOX_OK;
 }
 
+int feenox_problem_stiffness_add(element_t *e, unsigned int q, gsl_matrix *Kiq) {
 
-int feenox_problem_rhs_set(element_t *e, unsigned int q, double *value) {
+#ifdef HAVE_PETSC
+  double wdet = feenox_fem_compute_w_det_at_gauss(e, q);
+  gsl_matrix_scale(Kiq, wdet);
+  gsl_matrix_add(feenox.fem.Ki, Kiq);
+#endif
+
+  return FEENOX_OK;
+}
+
+int feenox_problem_rhs_add(element_t *e, unsigned int q, double *value) {
 
 #ifdef HAVE_PETSC
   for (unsigned int g = 0; g < feenox.pde.dofs; g++) {
     gsl_vector_set(feenox.fem.vec_f, g, value[g]);
   }  
 
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   gsl_matrix *H_Gc = feenox_fem_compute_H_Gc_at_gauss(e, q, feenox.pde.mesh->integration);
   feenox_call(feenox_blas_Atb_accum(H_Gc, feenox.fem.vec_f, wdet, feenox.fem.bi));
 #endif
 
   return FEENOX_OK;
 }
+

src/pdes/fem.c

@@ -318,7 +318,7 @@ inline material_t *feenox_fem_get_material(element_t *e) {
   return (e->physical_group != NULL) ? e->physical_group->material : NULL;
 }
 
-inline double feenox_fem_compute_w_det_at_gauss(element_t *e, unsigned int q, int integration) {
+inline double feenox_fem_compute_w_det_at_gauss_integration(element_t *e, unsigned int q, int integration) {
 
   if (feenox.fem.cache_J == 0 && e->type != feenox.fem.current_gauss_type) {
     feenox_fem_elemental_caches_reset();
@@ -609,7 +609,7 @@ inline gsl_matrix *feenox_fem_compute_H_Gc_at_gauss(element_t *e, unsigned int q
   return e->type->H_Gc[q];  
 }
 
-inline gsl_matrix *feenox_fem_compute_B_at_gauss(element_t *e, unsigned int q, int integration) {
+inline gsl_matrix *feenox_fem_compute_B_at_gauss_integration(element_t *e, unsigned int q, int integration) {
 
   if (feenox.fem.cache_J == 0 && e->type != feenox.fem.current_gauss_type) {
     feenox_fem_elemental_caches_reset();
@@ -638,7 +638,7 @@ inline gsl_matrix *feenox_fem_compute_B_G_at_gauss(element_t *e, unsigned int q,
   unsigned int G = feenox.pde.dofs;
   if (G == 1) {
     // for G=1 B_G = B
-    return feenox_fem_compute_B_at_gauss(e, q, integration);
+    return feenox_fem_compute_B_at_gauss_integration(e, q, integration);
   }
   gsl_matrix ***B_G = (feenox.fem.cache_B) ? &e->B_G : &feenox.fem.B_Gi;
   if ((*B_G) == NULL) {
@@ -654,7 +654,7 @@ inline gsl_matrix *feenox_fem_compute_B_G_at_gauss(element_t *e, unsigned int q,
     return (*B_G)[q];
   }
 
-  gsl_matrix *B = feenox_fem_compute_B_at_gauss(e, q, integration);
+  gsl_matrix *B = feenox_fem_compute_B_at_gauss_integration(e, q, integration);
   for (unsigned int d = 0; d < D; d++) {
     size_t Gd = G*d;
     for (unsigned int j = 0; j < J; j++) {

src/pdes/gradient.c

@@ -177,7 +177,7 @@ int feenox_problem_gradient_compute_at_element(element_t *e, mesh_t *mesh) {
       } else {
         gsl_matrix_set_zero(e->dphidx_gauss[q]);
       }
-      gsl_matrix *B = feenox_fem_compute_B_at_gauss(e, q, mesh->integration);
+      gsl_matrix *B = feenox_fem_compute_B_at_gauss_integration(e, q, mesh->integration);
 
       // aca habria que hacer una matriz con los phi globales
       // (de j y g, que de paso no depende de q asi que se podria hacer afuera del for de q)

src/pdes/laplace/bc.c

@@ -76,13 +76,13 @@ int feenox_problem_bc_set_laplace_derivative(bc_data_t *this, element_t *e, unsi
   // TODO: cache if neither space nor temperature dependent
   double *x = feenox_fem_compute_x_at_gauss_if_needed_and_update_var(e, q, feenox.pde.mesh->integration, this->space_dependent);
   double derivative = feenox_expression_eval(&this->expr);
-  feenox_call(feenox_problem_rhs_set(e, q, &derivative));
+  feenox_call(feenox_problem_rhs_add(e, q, &derivative));
 
   if (this->nonlinear) {
     double phi = feenox_function_eval(feenox.pde.solution[0], x);
     double dderivativedphi = feenox_expression_derivative_wrt_function(&this->expr, feenox.pde.solution[0], phi);
 
-    double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+    double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
     gsl_matrix *H_Gc = feenox_fem_compute_H_Gc_at_gauss(e, q, feenox.pde.mesh->integration);
     // mind the positive sign!
     feenox_call(feenox_blas_BtB_accum(H_Gc, +wdet*dderivativedphi, feenox.fem.Jbi));

src/pdes/laplace/bulk.c

@@ -26,20 +26,20 @@ int feenox_problem_build_volumetric_gauss_point_laplace(element_t *e, unsigned i
 
 #ifdef HAVE_PETSC
 
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
-  gsl_matrix *B = feenox_fem_compute_B_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss(e, q);
+  gsl_matrix *B = feenox_fem_compute_B_at_gauss(e, q);
 
   // laplace stiffness matrix Ki += wdet * Bt*B
   feenox_call(feenox_blas_BtB_accum(B, wdet, feenox.fem.Ki));
-
+  
   // the material is needed for either RHS and/or mass
   material_t *material = feenox_fem_get_material(e);
 
   // right-hand side
   double *x = feenox_fem_compute_x_at_gauss_if_needed(e, q, feenox.pde.mesh->integration, laplace.space_dependent_source || laplace.space_dependent_mass);
   if (laplace.f.defined) {
     double f = laplace.f.eval(&laplace.f, x, material);
-    feenox_call(feenox_problem_rhs_set(e, q, &f));
+    feenox_call(feenox_problem_rhs_add(e, q, &f));
   }
 
   if (feenox.pde.has_jacobian) {

src/pdes/mechanical/bc.c

@@ -339,7 +339,7 @@ int feenox_problem_bc_set_mechanical_normal_stress(bc_data_t *this, element_t *e
   }
 
   //feenox_call(feenox_problem_bc_natural_set(e, v, t));
-  feenox_call(feenox_problem_rhs_set(e, q, t));
+  feenox_call(feenox_problem_rhs_add(e, q, t));
 
 #endif
 
@@ -362,7 +362,7 @@ int feenox_problem_bc_set_mechanical_traction(bc_data_t *this, element_t *e, uns
 //  printf("%g\n", t[bc_data->dof]);
 //  feenox_call(feenox_problem_bc_natural_set(e, v, t));
 
-  feenox_call(feenox_problem_rhs_set(e, q, t));  
+  feenox_call(feenox_problem_rhs_add(e, q, t));  
 
 #endif
 
@@ -381,7 +381,7 @@ int feenox_problem_bc_set_mechanical_force(bc_data_t *this, element_t *e, unsign
     feenox_call(feenox_physical_group_compute_volume(e->physical_group, feenox.pde.mesh));
   }
   t[this->dof] = feenox_expression_eval(&this->expr) / e->physical_group->volume;
-  feenox_call(feenox_problem_rhs_set(e, q, t));  
+  feenox_call(feenox_problem_rhs_add(e, q, t));  
 
 #endif
 

src/pdes/mechanical/bulk.c

@@ -52,7 +52,7 @@ int feenox_problem_build_volumetric_gauss_point_mechanical(element_t *e, unsigne
     mechanical.compute_C(x, feenox_fem_get_material(e));
   }
 
-  gsl_matrix *dhdx = feenox_fem_compute_B_at_gauss(e, q, feenox.pde.mesh->integration);
+  gsl_matrix *dhdx = feenox_fem_compute_B_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   for (unsigned int j = 0; j < mechanical.n_nodes; j++) {
     // TODO: virtual methods? they cannot be inlined...
     if (mechanical.variant == variant_full) {
@@ -91,7 +91,7 @@ int feenox_problem_build_volumetric_gauss_point_mechanical(element_t *e, unsigne
   }
 
   // wdet
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
 
   // volumetric force densities
   if (mechanical.f_x.defined || mechanical.f_y.defined || mechanical.f_z.defined) {

src/pdes/modal/bulk.c

@@ -50,7 +50,7 @@ int feenox_problem_build_volumetric_gauss_point_modal(element_t *e, unsigned int
   }
 
   // TODO: unify with mechanical!
-  gsl_matrix *dhdx = feenox_fem_compute_B_at_gauss(e, q, feenox.pde.mesh->integration);
+  gsl_matrix *dhdx = feenox_fem_compute_B_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   for (unsigned int j = 0; j < modal.n_nodes; j++) {
     // TODO: virtual methods? they cannot be inlined...
     if (modal.variant == variant_full) {
@@ -89,7 +89,7 @@ int feenox_problem_build_volumetric_gauss_point_modal(element_t *e, unsigned int
   }
 
   // wdet
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
 
   // elemental stiffness B'*C*B
   feenox_call(feenox_blas_BtCB_accum(modal.B, modal.C, modal.CB, wdet, feenox.fem.Ki));

src/pdes/neutron_diffusion/bc.c

@@ -133,7 +133,7 @@ int feenox_problem_bc_set_neutron_diffusion_vacuum(bc_data_t *this, element_t *e
   double coeff = (this->expr.items != NULL) ? feenox_expression_eval(&this->expr) : 0.5;
 
   // TODO: convenience call like  feenox_problem_rhs_set()?
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   gsl_matrix *H = feenox_fem_compute_H_Gc_at_gauss(e, q, feenox.pde.mesh->integration);
   feenox_call(feenox_blas_BtB_accum(H, wdet*coeff, feenox.fem.Ki));
 
@@ -151,7 +151,7 @@ int feenox_problem_bc_set_neutron_diffusion_current(bc_data_t *this, element_t *
   feenox_fem_compute_x_at_gauss_if_needed_and_update_var(e, q, feenox.pde.mesh->integration, this->space_dependent);
   double *current = calloc(neutron_sn.groups, sizeof(double)); 
   current[this->dof] = feenox_expression_eval(&this->expr);
-  feenox_call(feenox_problem_rhs_set(e, q, current));
+  feenox_call(feenox_problem_rhs_add(e, q, current));
   feenox_free(current);
 
 #endif

src/pdes/neutron_diffusion/bulk.c

@@ -127,7 +127,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_diffusion(element_t *e,
   }
 
   // elemental stiffness for the diffusion term B'*D*B
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   gsl_matrix *B = feenox_fem_compute_B_G_at_gauss(e, q, feenox.pde.mesh->integration);
   feenox_call(feenox_blas_BtCB(B, neutron_diffusion.D_G, neutron_diffusion.DB, wdet, neutron_diffusion.Li));
 

src/pdes/neutron_sn/bulk.c

@@ -176,7 +176,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
   double tau = feenox_var_value(neutron_sn.sn_alpha) * e->type->size(e);
 
   gsl_matrix *H_G = feenox_fem_compute_H_Gc_at_gauss(e, q, feenox.pde.mesh->integration);
-  gsl_matrix *B = feenox_fem_compute_B_at_gauss(e, q, feenox.pde.mesh->integration);
+  gsl_matrix *B = feenox_fem_compute_B_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   feenox_call(gsl_matrix_memcpy(neutron_sn.P, H_G));
   for (unsigned int j = 0; j < neutron_sn.n_nodes; j++) {
     int MGj = MG*j;
@@ -192,7 +192,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
   }
 
   // petrov-stabilized leakage term
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   gsl_matrix *B_G = feenox_fem_compute_B_G_at_gauss(e, q, feenox.pde.mesh->integration);
   feenox_call(feenox_blas_PtCB(neutron_sn.P, neutron_sn.D, B_G, neutron_sn.DB, wdet, neutron_sn.Li));
 

src/pdes/neutron_sn/post.c

@@ -47,7 +47,7 @@ int feenox_problem_solve_post_neutron_sn(void) {
         // for now I'd rather perform the computation myself instead of 
         // calling feenox_mesh_integral_over_element() because that wastes a lot of fem calls
         for (unsigned int q = 0; q < element->type->gauss[integration].Q; q++) {
-          double wdet = feenox_fem_compute_w_det_at_gauss(element, q, integration);
+          double wdet = feenox_fem_compute_w_det_at_gauss_integration(element, q, integration);
           for (unsigned int j = 0; j < element->type->nodes; j++) {
             double h = gsl_matrix_get(element->type->gauss[integration].H_c[q], 0, j);
             for (unsigned int g = 0; g < feenox.pde.dofs; g++) {

src/pdes/thermal/bc.c

@@ -105,14 +105,14 @@ int feenox_problem_bc_set_thermal_heatflux(bc_data_t *this, element_t *e, unsign
   // TODO: cache if neither space nor temperature dependent
   double *x = feenox_fem_compute_x_at_gauss_if_needed_and_update_var(e, q, feenox.pde.mesh->integration, this->space_dependent);
   double power = feenox_expression_eval(&this->expr);
-  feenox_call(feenox_problem_rhs_set(e, q, &power));
+  feenox_call(feenox_problem_rhs_add(e, q, &power));
 
   if (this->nonlinear) {
     double T = feenox_function_eval(feenox.pde.solution[0], x);
     double dqdT = feenox_expression_derivative_wrt_function(&this->expr, feenox.pde.solution[0], T);
     // mind the positive sign!
     gsl_matrix *H = feenox_fem_compute_H_Gc_at_gauss(e, q, feenox.pde.mesh->integration);
-    double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+    double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
     feenox_call(feenox_blas_BtB_accum(H, +wdet*dqdT, feenox.fem.Jbi));
   }
 
@@ -149,11 +149,11 @@ int feenox_problem_bc_set_thermal_convection(bc_data_t *this, element_t *e, unsi
 
   // the h*Tref goes to b
   double rhs = h*Tref;
-  feenox_call(feenox_problem_rhs_set(e, q, &rhs));
+  feenox_call(feenox_problem_rhs_add(e, q, &rhs));
 
   // TODO: the h*T goes directly to the stiffness matrix
   // this is not efficient because if h depends on t or T we might need to re-build the whole K
-  double wdet = feenox_fem_compute_w_det_at_gauss(e, q, feenox.pde.mesh->integration);
+  double wdet = feenox_fem_compute_w_det_at_gauss_integration(e, q, feenox.pde.mesh->integration);
   gsl_matrix *H = feenox_fem_compute_H_Gc_at_gauss(e, q, feenox.pde.mesh->integration);
   feenox_call(feenox_blas_BtB_accum(H, +wdet*h, feenox.fem.Ki));