Force BCs (#27)

* Initial commit

* Debug of lattice issue

* FIxed error by mapping boundary attributes

* Added yaml file

* Working example

* Changed default so that solver passes tests

* Works

* Removed unnecessary comments

* Cleaned up the function handling

* Removed unnecessary file outputs

* Making sure that checks are running

* ROM boundary tweak

* removed print statement

* Removed type_idx and loop

* Fixed double indexing

* Cleaner version

* Changed requested changes

examples/linelast/CMakeLists.txt

@@ -6,6 +6,7 @@ file(COPY linelast.simpleL.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast
 file(COPY config/linelast.simpleL.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/config)
 file(COPY linelast.lattice.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast)
 file(COPY config/linelast.lattice.h5 DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/config)
+file(COPY linelast.force_cantilever.yml DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast)
 
 file(COPY meshes/beam-tri.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)
 file(COPY meshes/beam-tet.mesh DESTINATION ${CMAKE_BINARY_DIR}/examples/linelast/meshes)

examples/linelast/linelast.force_cantilever.yml

@@ -0,0 +1,74 @@
+main:
+  mode: single_run
+  use_rom: true
+  solver: linelast
+
+mesh:
+  type: component-wise
+  component-wise:
+    global_config: "config/linelast.lattice.h5"
+    components:
+      - name: "joint2D"
+        file: "meshes/joint2D.mesh"
+      - name: "rod2D_H"
+        file: "meshes/rod2D_H.mesh"
+      - name: "rod2D_V"
+        file: "meshes/rod2D_V.mesh"
+
+domain-decomposition:
+  type: interior_penalty
+
+discretization:
+  order: 1
+  full-discrete-galerkin: true
+
+visualization:
+  enabled: false
+  unified_paraview: true
+  file_path:
+    prefix: paraview_output
+
+parameterized_problem:
+  name: linelast_force_cantilever
+
+single_run:
+  linelast_force_cantilever:
+    rforce_f: 1.0
+
+sample_generation:
+  maximum_number_of_snapshots: 400
+  component_sampling: false
+  file_path:
+    prefix: "linelast_force_cantilever"
+  parameters:
+    - key: single_run/linelast_force/rforce_f
+      type: double
+      sample_size: 1
+      minimum: 1.0
+      maximum: 1.0
+
+basis:
+  prefix: "linelast_force_cantilever"
+  tags:
+    - name: "joint2D"
+      number_of_basis: 12
+    - name: "rod2D_H"
+      number_of_basis: 8
+    - name: "rod2D_V"
+      number_of_basis: 9
+  svd:
+    save_spectrum: true
+    update_right_sv: false
+  visualization:
+    enabled: false
+
+model_reduction:
+  rom_handler_type: mfem
+  subdomain_training: universal
+  linear_solver_type: direct
+  save_operator:
+    level: component
+    prefix: "linelast.rom_elem"
+  compare_solution:
+    enabled: true
+

include/linelast_solver.hpp

@@ -52,6 +52,9 @@ class LinElastSolver : public MultiBlockSolver
    // Initial positions
    VectorFunctionCoefficient *init_x = NULL;
 
+   // Boundary condition types
+   Array<LinElastProblem::BoundaryType> type_idx;
+
 public:
    LinElastSolver();
 

include/parameterized_problem.hpp

@@ -85,6 +85,13 @@ extern double rdisp_f;
 void init_disp(const Vector &x, Vector &u);
 }
 
+namespace linelast_force
+{
+extern double rforce_f;
+
+void tip_force(const Vector &x, Vector &u);
+}
+
 }
 
 class ParameterizedProblem
@@ -240,6 +247,12 @@ class LinElastDispLattice : public LinElastProblem
    LinElastDispLattice();
 };
 
+class LinElastForceCantilever : public LinElastProblem
+{
+public:
+   LinElastForceCantilever();
+};
+
 ParameterizedProblem* InitParameterizedProblem();
 
 #endif

src/linelast_solver.cpp

@@ -75,6 +75,9 @@ void LinElastSolver::SetupBCVariables()
    bdr_coeffs.SetSize(numBdr);
    bdr_coeffs = NULL;
 
+   type_idx.SetSize(numBdr);
+   type_idx = LinElastProblem::BoundaryType::DIRICHLET;
+
    lambda_c.SetSize(numSub);
    lambda_c = NULL;
 
@@ -145,7 +148,10 @@ bool LinElastSolver::BCExistsOnBdr(const int &global_battr_idx)
 {
    assert((global_battr_idx >= 0) && (global_battr_idx < global_bdr_attributes.Size()));
    assert(bdr_coeffs.Size() == global_bdr_attributes.Size());
-   return (bdr_coeffs[global_battr_idx]);
+   if(type_idx[global_battr_idx] == LinElastProblem::BoundaryType::NEUMANN)
+      return false;
+   else
+      return (bdr_coeffs[global_battr_idx]);
 }
 
 void LinElastSolver::BuildRHSOperators()
@@ -172,10 +178,19 @@ void LinElastSolver::SetupRHSBCOperators()
          int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]);
          if (idx < 0)
             continue;
+
+         if (type_idx[b] == LinElastProblem::BoundaryType::NEUMANN)
+         {
+            bs[m]->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(*bdr_coeffs[b]), *bdr_markers[b]);
+         }  
+
          if (!BCExistsOnBdr(b))
             continue;
 
-         bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*bdr_coeffs[b], *lambda_c[m], *mu_c[m], alpha, kappa), *bdr_markers[b]);
+         if (type_idx[b] == LinElastProblem::BoundaryType::DIRICHLET) // Default behavior
+         { 
+            bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*bdr_coeffs[b], *lambda_c[m], *mu_c[m], alpha, kappa), *bdr_markers[b]);
+         }
       }
    }
 }
@@ -438,21 +453,19 @@ void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem)
       mu_c[i] = new ConstantCoefficient(mu_i);
    }
 
-   // Set BCs
-   for (int b = 0; b < problem->battr.Size(); b++)
+   // Set BCs, the switch on BC type is done inside SetupRHSBCOperators
+   for (int b = 0; b < global_bdr_attributes.Size(); b++)
    {
-      switch (problem->bdr_type[b])
-      {
-      case LinElastProblem::BoundaryType::NEUMANN: break;
-      case LinElastProblem::BoundaryType::ZERO: break;
+      int ti = problem->battr.Find(global_bdr_attributes[b]);
 
-      default:
-      case LinElastProblem::BoundaryType::DIRICHLET:
+      if (ti >=0)
       {
-         assert(problem->vector_bdr_ptr[b]);
-         AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]);
-         break;
-      }
+         type_idx[b] = (LinElastProblem::BoundaryType)problem->bdr_type[ti];
+         if (problem->bdr_type[ti] == LinElastProblem::BoundaryType::DIRICHLET || problem->bdr_type[ti] == LinElastProblem::BoundaryType::NEUMANN)
+         {
+            assert(problem->vector_bdr_ptr[ti]);
+            AddBCFunction(*(problem->vector_bdr_ptr[ti]), problem->battr[ti]);
+         }
       }
    }
 
@@ -525,8 +538,8 @@ void LinElastSolver::BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp)
          bdr_marker = 0;
          bdr_marker[comp->bdr_attributes[b] - 1] = 1;
          BilinearForm a_comp(fes_comp[c]);
-
          a_comp.AddBdrFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]), alpha, kappa), bdr_marker);
+
          a_comp.Assemble();
          a_comp.Finalize();
 

src/parameterized_problem.cpp

@@ -187,6 +187,18 @@ void init_disp_lcantilever(const Vector &x, Vector &u)
 
 }  // namespace linelast_disp
 
+namespace linelast_force
+{
+
+double rforce_f;
+
+void tip_force(const Vector &x, Vector &f)
+{
+   f = 0.0;
+   f(f.Size()-1) = -1.0e-2* rforce_f;
+}
+
+}  // namespace linelast_force
 
 }  // namespace function_factory
 
@@ -287,6 +299,10 @@ ParameterizedProblem* InitParameterizedProblem()
    {
       problem = new LinElastDispLattice();
    }
+   else if (problem_name == "linelast_force_cantilever")
+   {
+      problem = new LinElastForceCantilever();
+   }
    else
    {
       mfem_error("Unknown parameterized problem name!\n");
@@ -577,15 +593,19 @@ LinElastDisp::LinElastDisp()
     : LinElastProblem()
 {
    // pointer to static function.
-   bdr_type.SetSize(2);
-   battr.SetSize(2);
-   vector_bdr_ptr.SetSize(2);
+   bdr_type.SetSize(3);
+   battr.SetSize(3);
+   vector_bdr_ptr.SetSize(3);
    for (size_t i = 0; i < vector_bdr_ptr.Size(); i++)
    {
    bdr_type[i] = LinElastProblem::DIRICHLET;
    battr[i] = i+1;
    vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp);
    }
+
+   battr[2] = 3;
+   bdr_type[2] = LinElastProblem::ZERO;
+   vector_bdr_ptr[2] = NULL;
 
    // Set materials
    general_scalar_ptr.SetSize(2);
@@ -615,15 +635,19 @@ LinElastDispLCantilever::LinElastDispLCantilever()
     : LinElastProblem()
 {
    // pointer to static function.
-   bdr_type.SetSize(2);
-   battr.SetSize(2);
-   vector_bdr_ptr.SetSize(2);
-   for (size_t i = 0; i < vector_bdr_ptr.Size(); i++)
+   bdr_type.SetSize(3);
+   battr.SetSize(3);
+   vector_bdr_ptr.SetSize(3);
+   for (size_t i = 0; i < 2; i++)
    {
-   bdr_type[i] = LinElastProblem::DIRICHLET;
-   battr[i] = i+1;
-   vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp_lcantilever);
+      battr[i] = i+1;
+      bdr_type[i] = LinElastProblem::DIRICHLET;
+      vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp_lcantilever);
    }
+
+   battr[2] = 3;
+   bdr_type[2] = LinElastProblem::ZERO;
+   vector_bdr_ptr[2] = NULL;
 
    // Set materials
    general_scalar_ptr.SetSize(2);
@@ -652,16 +676,20 @@ LinElastDispLattice::LinElastDispLattice()
     : LinElastProblem()
 {
    // pointer to static function.
-   bdr_type.SetSize(2);
-   battr.SetSize(2);
-   vector_bdr_ptr.SetSize(2);
-   for (size_t i = 0; i < vector_bdr_ptr.Size(); i++)
-   {
-   bdr_type[i] = LinElastProblem::DIRICHLET;
-   battr[i] = i+1;
-   vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp);
+   bdr_type.SetSize(3);
+   battr.SetSize(3);
+   vector_bdr_ptr.SetSize(3);
+  for (size_t i = 0; i < 2; i++)
+      {
+      battr[i] = i+1;
+      bdr_type[i] = LinElastProblem::DIRICHLET;
+      vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp);
    }
-
+
+   battr[2] = 3;
+   bdr_type[2] = LinElastProblem::ZERO;
+   vector_bdr_ptr[2] = NULL;
+
    // Set materials
    general_scalar_ptr.SetSize(2);
    general_scalar_ptr[0] = function_factory::linelast_problem::lambda;
@@ -685,3 +713,49 @@ LinElastDispLattice::LinElastDispLattice()
    general_vector_ptr[0] = NULL;
 
 }
+
+
+LinElastForceCantilever::LinElastForceCantilever()
+    : LinElastProblem()
+{
+   // pointer to static function.
+   bdr_type.SetSize(3);
+   battr.SetSize(3);
+   vector_bdr_ptr.SetSize(3);
+
+   // Fixed end of cantilever
+   battr[0] = 1;
+   bdr_type[0] = LinElastProblem::DIRICHLET;
+   vector_bdr_ptr[0] = &(function_factory::linelast_disp::init_disp);
+
+   // Free end of cantilever
+   battr[1] = 2;
+   bdr_type[1] = LinElastProblem::NEUMANN;
+   vector_bdr_ptr[1] = &(function_factory::linelast_force::tip_force);
+
+   battr[2] = 3;
+   bdr_type[2] = LinElastProblem::ZERO;
+   vector_bdr_ptr[2] = NULL;
+
+   // Set materials
+   general_scalar_ptr.SetSize(2);
+   general_scalar_ptr[0] = function_factory::linelast_problem::lambda;
+   general_scalar_ptr[1] = function_factory::linelast_problem::mu;
+
+   // Default values.
+   function_factory::linelast_force::rforce_f = 1.0;
+   function_factory::linelast_problem::_lambda = 1.0;
+   function_factory::linelast_problem::_mu = 1.0;
+
+   param_map["rforce_f"] = 0;
+   param_map["lambda"] = 1;
+   param_map["mu"] = 2;
+
+   param_ptr.SetSize(3);
+   param_ptr[0] = &(function_factory::linelast_force::rforce_f);
+   param_ptr[1] = &(function_factory::linelast_problem::_lambda);
+   param_ptr[2] = &(function_factory::linelast_problem::_mu);
+
+   general_vector_ptr.SetSize(1);
+   general_vector_ptr[0] = NULL;
+}