nomenclature and laplace rhs like H^T*vec_f

src/mesh/fem.c

@@ -582,23 +582,27 @@ int feenox_mesh_compute_x_at_gauss(element_t *e, unsigned int q, int integration
 
 int feenox_mesh_compute_H_at_gauss(element_t *e, unsigned int q, int integration) {
   // TODO: esta H no depende del elemento, depende solamente de la cantidad de grados de libertad del problema
+  //       y de la cantidad de nodos entonces habria que mandarla a feenox.pde en lugar de hacerla por elemento
   if (e->H == NULL) {
     feenox_check_alloc(e->H = calloc(e->type->gauss[integration].Q, sizeof(gsl_matrix *)));
   }
 
-  unsigned int dofs = feenox.pde.dofs;
+  unsigned int G = feenox.pde.dofs;
 
   if (e->H[q] == NULL) {
-    feenox_check_alloc(e->H[q] = gsl_matrix_calloc(dofs, dofs*e->type->nodes));
+    feenox_check_alloc(e->H[q] = gsl_matrix_calloc(G, G*e->type->nodes));
   } else {
     return FEENOX_OK;
   }
 
-  for (unsigned int g = 0; g < dofs; g++) {
+  for (unsigned int g = 0; g < G; g++) {
     for (unsigned int j = 0; j < e->type->nodes; j++) {
-      gsl_matrix_set(e->H[q], g, dofs*j+g, e->type->gauss[integration].h[q][j]);
+      gsl_matrix_set(e->H[q], g, G*j+g, e->type->gauss[integration].h[q][j]);
     }
   }
+
+//  feenox_debug_print_gsl_matrix(e->H[q], stdout);
+
 
   return FEENOX_OK;  
 }
@@ -624,16 +628,18 @@ int feenox_mesh_compute_B_at_gauss(element_t *e, unsigned int q, int integration
   // mondiola es una transpose! es la misma pero venimos transpuestos nosotros!
   // TODO: if dofs == 1 apuntar a e->dhdx y ya
   for (unsigned int d = 0; d < e->type->dim; d++) {
-    size_t dofs_times_d = feenox.pde.dofs*d;
+    size_t Gd = feenox.pde.dofs*d;
     for (unsigned int j = 0; j < e->type->nodes; j++) {
-      size_t dofs_times_j = feenox.pde.dofs*j;
-      double dhdx = gsl_matrix_get(e->dhdx[q], j, d);
+      size_t Gj = feenox.pde.dofs*j;
+      double dhjdx = gsl_matrix_get(e->dhdx[q], j, d);
       for (unsigned int g = 0; g < feenox.pde.dofs; g++) {
-        gsl_matrix_set(e->B[q], dofs_times_d+g, dofs_times_j+g, dhdx);
+        gsl_matrix_set(e->B[q], Gd+g, Gj+g, dhjdx);
       }
     }
   }
-
+
+//  feenox_debug_print_gsl_matrix(e->B[q], stdout);
+
   return FEENOX_OK;
 }
 

src/pdes/laplace/bc.c

@@ -69,22 +69,22 @@ int feenox_problem_bc_set_laplace_phi(bc_data_t *this, element_t *e, size_t j_gl
   return FEENOX_OK;
 }
 
-int feenox_problem_bc_set_laplace_derivative(bc_data_t *this, element_t *e, unsigned int v) {
+int feenox_problem_bc_set_laplace_derivative(bc_data_t *this, element_t *e, unsigned int q) {
 
 #ifdef HAVE_PETSC
 
   // TODO: cache if neither space nor temperature dependent
-  double *x = feenox_problem_bc_natural_x(e, this, v);
+  double *x = feenox_problem_bc_natural_x(e, this, q);
   double derivative = feenox_expression_eval(&this->expr);
-  feenox_call(feenox_problem_bc_natural_set(e, v, &derivative));
+  feenox_call(feenox_problem_bc_natural_set(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);
     // TODO: axisymmetric
-    double w = e->w[v];
+    double w = e->w[q];
     // mind the positive sign!
-    feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, +w*dderivativedphi, e->H[v], e->H[v], 1.0, feenox.pde.Jbi));
+    feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, +w*dderivativedphi, e->H[q], e->H[q], 1.0, feenox.pde.Jbi));
   }
 
 #endif

src/pdes/laplace/bulk.c

@@ -1,7 +1,7 @@
 /*------------ -------------- -------- --- ----- ---   --       -            -
  *  feenox's routines for Laplace's equation: bulk elements
  *
- *  Copyright (C) 2021--2022 jeremy theler
+ *  Copyright (C) 2021--2023 jeremy theler
  *
  *  This file is part of FeenoX <https://www.seamplex.com/feenox>.
  *
@@ -22,35 +22,28 @@
 #include "feenox.h"
 #include "laplace.h"
 
-int feenox_problem_build_volumetric_gauss_point_laplace(element_t *e, unsigned int v) {
+int feenox_problem_build_volumetric_gauss_point_laplace(element_t *e, unsigned int q) {
 
 #ifdef HAVE_PETSC
 
-  feenox_call(feenox_mesh_compute_wHB_at_gauss(e, v));
-  double *x = feenox_mesh_compute_x_if_needed(e, v, laplace.space_dependent_source || laplace.space_dependent_mass);
+  feenox_call(feenox_mesh_compute_wHB_at_gauss(e, q));
+  double *x = feenox_mesh_compute_x_if_needed(e, q, laplace.space_dependent_source || laplace.space_dependent_mass);
 
   // TODO: axisymmetric
 //  r_for_axisymmetric = feenox_compute_r_for_axisymmetric(this, v);
-  double r_for_axisymmetric = 1;
-  double w = e->w[v] * r_for_axisymmetric;
+//  double r_for_axisymmetric = 1;
+  double w = e->w[q];
 
   // laplace stiffness matrix Bt*B
   // note we don't allow any coefficient in the laplacian term
-  feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, w, e->B[v], e->B[v], 1.0, feenox.pde.Ki));
+  feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, w, e->B[q], e->B[q], 1.0, feenox.pde.Ki));
 
-  material_t *material = NULL;
-  if (e->physical_group != NULL && e->physical_group->material != NULL) {
-    material = e->physical_group->material;
-  }
+  material_t *material = feenox_mesh_get_material(e);
 
   // right-hand side
-  // TODO: there should be a way to use BLAS instead of a for loop
   if (laplace.f.defined) {
-    double f = laplace.f.eval(&laplace.f, x, material);
-    unsigned int j = 0;
-    for (j = 0; j < e->type->nodes; j++) {
-      gsl_vector_add_to_element(feenox.pde.bi, j, w * e->type->gauss[feenox.pde.mesh->integration].h[v][j] * f);
-    }
+    gsl_vector_set(laplace.vec_f, 0, laplace.f.eval(&laplace.f, x, material));
+    feenox_call(gsl_blas_dgemv(CblasTrans, w, e->H[q], laplace.vec_f, 1.0, feenox.pde.bi));
   }
 
   if (feenox.pde.has_jacobian) {
@@ -67,7 +60,7 @@ int feenox_problem_build_volumetric_gauss_point_laplace(element_t *e, unsigned i
       feenox_push_error_message("no alpha found");
       return FEENOX_ERROR;
     }
-    feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, w * alpha, e->H[v], e->H[v], 1.0, feenox.pde.Mi));
+    feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, w * alpha, e->H[q], e->H[q], 1.0, feenox.pde.Mi));
   }
 
 

src/pdes/laplace/init.c

@@ -96,6 +96,10 @@ int feenox_problem_init_runtime_laplace(void) {
     }
   }
 
+  // allocate vector for computing the RHS as a matrix-vector product H^T*vec_f
+  feenox_check_alloc(laplace.vec_f = gsl_vector_alloc(1));
+
+
   laplace.alpha.uniform = feenox_expression_depends_on_space(laplace.alpha.dependency_variables);
   laplace.alpha.non_linear      = feenox_expression_depends_on_function(laplace.alpha.dependency_functions,  feenox.pde.solution[0]);
 

src/pdes/laplace/laplace.h

@@ -39,6 +39,9 @@ struct laplace_t {
   int space_dependent_bc;
   int phi_dependent_bc;
 
+  // right-hand-side vector of size one for the product H^T*f
+  gsl_vector *vec_f;
+
   // caches for uniform properties
   struct {
     double f;

src/pdes/laplace/methods.h

@@ -1,7 +1,7 @@
 /*------------ -------------- -------- --- ----- ---   --       -            -
  *  feenox's routines for Laplace's equation: virtual methods
  *
- *  Copyright (C) 2021 jeremy theler
+ *  Copyright (C) 2022-2023 jeremy theler
  *
  *  This file is part of FeenoX <https://www.seamplex.com/feenox>.
  *

src/pdes/neutron_diffusion/bulk.c

@@ -141,10 +141,10 @@ int feenox_problem_build_volumetric_gauss_point_neutron_diffusion(element_t *e,
   }
 
   // for source-driven problems
-  //   K = Ki + Ai - Xi
+  //   Ki = Li + Ai - Xi
   // for criticallity problems
-  //   K = Ki + Ai
-  //   M = Xi
+  //   Ki = Li + Ai
+  //   Mi = Xi
   gsl_matrix_add(neutron_diffusion.Li, neutron_diffusion.Ai);
   if (neutron_diffusion.has_fission) {
     if (neutron_diffusion.has_sources) {

tests/laplace2d.sh

@@ -13,5 +13,9 @@ checkpde laplace
 checkgmsh
 
 gmsh -2 ${dir}/square-centered.geo || exit $?
+
 answer laplace-square.fee "-1 -2.68735 -1"
 exitifwrong $?
+
+answer poisson-square.fee "0.2949 0.2044 -0.2044"
+exitifwrong $?

tests/poisson-square.fee

@@ -0,0 +1,7 @@
+PROBLEM laplace 2d
+READ_MESH square-centered.msh
+BC zero  phi=0 PHYSICAL_GROUP left PHYSICAL_GROUP right PHYSICAL_GROUP bottom PHYSICAL_GROUP up
+f = 1
+SOLVE_PROBLEM
+WRITE_RESULTS
+PRINT %.4f phi(0,0) dphidx(-0.5,-0.5) dphidy(0.5,0.5) SEP " "