GismoLinearElasticityModel.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <GismoPdeModel.hpp>
 
using namespace std::string_literals;
 
namespace iganet {
 
namespace webapp {
 
template <short_t d, typename T>
class GismoLinearElasticityModel : public GismoPdeModel<d, T> {
 
  static_assert(d >= 1 && d <= 3, "Spatial dimension must be between 1 and 3");
 
private:
  using Base = GismoPdeModel<d, T>;
 
  gismo::gsMultiBasis<T> basis_;
 
  gismo::gsBoundaryConditions<T> bc_;
 
  gismo::gsFunctionExpr<T> rhsFunc_;
  gismo::gsFunctionExpr<T> loadFunc_;
 
  std::array<gismo::gsFunctionExpr<T>, 2 * d> bcFunc_;
 
  std::array<gismo::condition_type::type, 2 * d> bcType_;
 
  T YoungsModulus_;
 
  T PoissonsRatio_;
 
  void solve() {
 
    // Setup assembler
    gismo::gsElasticityAssembler<T> assembler(Base::geo_, basis_, bc_, rhsFunc_);
    assembler.options().setReal("YoungsModulus", YoungsModulus_);
    assembler.options().setReal("PoissonsRatio", PoissonsRatio_);
    assembler.options().setInt("MaterialLaw", gismo::material_law::hooke);
    assembler.options().setInt("DirichletStrategy", gismo::dirichlet::elimination);
 
    // Initialize assembler
    assembler.assemble();
 
    // Solve system
    typename gismo::gsSparseSolver<T>::CGDiagonal solver(assembler.matrix());
    gismo::gsMatrix<T> solution(solver.solve(assembler.rhs()));
 
    // Extract solution
    assembler.constructSolution(solution, assembler.allFixedDofs(),
                                Base::solution_);
  }
 
public:
  GismoLinearElasticityModel() = delete;
 
  GismoLinearElasticityModel(const std::array<short_t, d> degrees,
                             const std::array<int64_t, d> ncoeffs,
                             const std::array<int64_t, d> npatches,
                             const std::array<T, d> dimensions)
      : Base(degrees, ncoeffs, npatches, dimensions), basis_(Base::geo_, true),
        YoungsModulus_(210e9), PoissonsRatio_(0.3), rhsFunc_("0", "0", "0", 3),
        loadFunc_("0", "0", "-1e5", 3) {
 
    // Set boundary conditions type and expression
    for (const auto &side : GismoBoundarySides<d>) {
      bcType_[side - 1] = gismo::condition_type::unknownType;
      bcFunc_[side - 1] = gismo::give(gismo::gsFunctionExpr<T>("0", "0", "0", 3));
      // bc_.addCondition(0, side, bcType_[side - 1], &bcFunc_[side - 1]);
    }
 
    bc_.addCondition(0, gismo::boundary::west, gismo::condition_type::dirichlet,
                     nullptr, 0);
    bc_.addCondition(0, gismo::boundary::west, gismo::condition_type::dirichlet,
                     nullptr, 1);
    bc_.addCondition(0, gismo::boundary::west, gismo::condition_type::dirichlet,
                     nullptr, 2);
 
    bc_.addCondition(0, gismo::boundary::east, gismo::condition_type::neumann,
                     &loadFunc_);
 
    // Set geometry
    bc_.setGeoMap(Base::geo_);
 
    // Regenerate solution
    solve();
  }
 
  ~GismoLinearElasticityModel() {}
 
  std::string getName() const override {
    return "GismoLinearElasticity" + std::to_string(d) + "d";
  }
 
  std::string getDescription() const override {
    return "G+Smo linear elasticity model in " + std::to_string(d) +
           " dimensions";
  };
 
  nlohmann::json getOptions() const override {
 
    nlohmann::json json = Base::getOptions();
 
    return json;
  }
 
  nlohmann::json getOutputs() const override {
    auto json = R"([{
           "name" : "Displacement",
           "description" : "Displament magnitude",
           "type" : 1},{
           "name" : "Displacement_x",
           "description" : "Displacement x-component",
           "type" : 1},{
           "name" : "Displacement_y",
           "description" : "Displacement x-component",
           "type" : 1},{
           "name" : "Displacement_z",
           "description" : "Displacement z-component",
           "type" : 1}])"_json;
 
    for (auto const &output : Base::getOutputs())
      json.push_back(output);
 
    return json;
  }
 
  nlohmann::json getParameters() const override {
 
    auto json = nlohmann::json::array();
    int uiid = 0;
 
    // Lambda expression to add a JSON entry
    auto add_json = [&json, &uiid]<typename Type, typename Value>(
                        const std::string &name, const std::string &label,
                        const std::string &group,
                        const std::string &description, const Type &type,
                        const Value &value) {
      nlohmann::json item;
      item["name"] = name;
      item["label"] = label;
      item["description"] = description;
      item["group"] = group;
      item["type"] = type;
      item["value"] = value;
      item["default"] = value;
      item["uuid"] = uiid++;
      json.push_back(item);
    };
 
    // Lambda expression to add a JSON entry with different default type
    auto add_json_default =
        [&json, &uiid]<typename Type, typename Value, typename DefaultValue>(
            const std::string &name, const std::string &label,
            const std::string &group, const std::string &description,
            const Type &type, const Value &value,
            const DefaultValue &defaultValue) {
          nlohmann::json item;
          item["name"] = name;
          item["label"] = label;
          item["description"] = description;
          item["group"] = group;
          item["type"] = type;
          item["value"] = value;
          item["default"] = defaultValue;
          item["uuid"] = uiid++;
          json.push_back(item);
        };
 
    add_json("YoungModulus", "Young", "", "Young's modulus", "float",
             YoungsModulus_);
    add_json("PoissonRatio", "Poisson", "", "Poisson's ratio", "float",
             PoissonsRatio_);
 
    //    add_json("rhs", "Right-hand side function",
    //             std::vector<std::string>{"text", "text", "text"},
    //             std::vector<std::string>{rhsFunc_.expression(0),
    //                                      rhsFunc_.expression(1),
    //                                      rhsFunc_.expression(2)});
 
    return json;
  }
 
  nlohmann::json updateAttribute(const std::string &patch,
                                 const std::string &component,
                                 const std::string &attribute,
                                 const nlohmann::json &json) override {
 
    nlohmann::json result = R"({})"_json;
 
    if (attribute == "YoungModulus") {
      if (!json.contains("data"))
        throw InvalidModelAttributeException();
      if (!json["data"].contains("YoungModulus"))
        throw InvalidModelAttributeException();
 
      YoungsModulus_ = json["data"]["YoungModulus"].get<T>();
    }
 
    else if (attribute == "PoissonRatio") {
      if (!json.contains("data"))
        throw InvalidModelAttributeException();
      if (!json["data"].contains("PoissonRatio"))
        throw InvalidModelAttributeException();
 
      PoissonsRatio_ = json["data"]["PoissonRatio"].get<T>();
    }
 
    else
      result = Base::updateAttribute(patch, component, attribute, json);
 
    // Solve updated problem
    solve();
 
    return result;
  }
 
  nlohmann::json eval(const std::string &patch, const std::string &component,
                      const nlohmann::json &json) const override {
 
    int patchIndex(-1);
 
    try {
      patchIndex = stoi(patch);
    } catch (...) {
      // Invalid patchIndex
      return R"({ INVALID REQUEST })"_json;
    }
 
    if (component == "Displacement" ||
        component == "Displacement_x" ||
        component == "Displacement_y" ||
        component == "Displacement_z") {
      
      nlohmann::json result;
 
      // degrees
      result["degrees"] = nlohmann::json::array();
      
      for (std::size_t i = 0; i < Base::solution_.patch(patchIndex).parDim(); ++i)
        result["degrees"].push_back(Base::solution_.patch(patchIndex).degree(i));
 
      // ncoeffs
      result["ncoeffs"] = nlohmann::json::array();
      
      if (auto bspline =
          dynamic_cast<const gismo::gsBSpline<T> *>(&Base::solution_.patch(patchIndex)))
        for (std::size_t i = 0; i < bspline->parDim(); ++i)
          result["ncoeffs"].push_back(bspline->basis().size(i));
      else if (auto bspline = dynamic_cast<const gismo::gsTensorBSpline<d, T> *>(
                                                                                 &Base::solution_.patch(patchIndex)))
        for (std::size_t i = 0; i < bspline->parDim(); ++i)
          result["ncoeffs"].push_back(bspline->basis().size(i));
      else
        return R"({ INVALID REQUEST })"_json;
      
      // nknots
      result["nknots"] = nlohmann::json::array();
 
      if (auto bspline =
          dynamic_cast<const gismo::gsBSpline<T> *>(&Base::solution_.patch(patchIndex)))
        for (std::size_t i = 0; i < bspline->parDim(); ++i)
          result["nknots"].push_back(bspline->knots(i).size());
      else if (auto bspline = dynamic_cast<const gismo::gsTensorBSpline<d, T> *>(
                                                                                 &Base::solution_.patch(patchIndex)))
        for (std::size_t i = 0; i < bspline->parDim(); ++i)
          result["nknots"].push_back(bspline->knots(i).size());
      else
        return R"({ INVALID REQUEST })"_json;
 
      // knots
      result["knots"] = nlohmann::json::array();
 
      if (auto bspline =
          dynamic_cast<const gismo::gsBSpline<T> *>(&Base::solution_.patch(patchIndex)))
        for (std::size_t i = 0; i < bspline->parDim(); ++i)
          result["knots"].push_back(bspline->knots(i));
      else if (auto bspline = dynamic_cast<const gismo::gsTensorBSpline<d, T> *>(
                                                                                 &Base::solution_.patch(patchIndex)))
        for (std::size_t i = 0; i < bspline->parDim(); ++i)
          result["knots"].push_back(bspline->knots(i));
      else
        return R"({ INVALID REQUEST })"_json;
      
      // coeffs
      gismo::gsMatrix<T> coeffs;
      if (component == "Displacement") {
        coeffs = Base::solution_.patch(patchIndex).coefs().rowwise().norm();        
      } else if (component == "Displacement_x") {
        coeffs = Base::solution_.patch(patchIndex).coefs().col(0);
      } else if (component == "Displacement_y") {
        coeffs = Base::solution_.patch(patchIndex).coefs().col(1);
      } else if (component == "Displacement_z") {
        coeffs = Base::solution_.patch(patchIndex).coefs().col(2);
      }
      result["coeffs"] = utils::to_json(coeffs, true, false);
      
      return result;      
    }
    
    else
      return Base::eval(patch, component, json);
  }
 
  void elevate(const nlohmann::json &json = NULL) override {
 
    bool geometry = true;
 
    if (json.contains("data"))
      if (json["data"].contains("num"))
        geometry = json["data"]["geometry"].get<bool>();
 
    if (geometry) {
      // Elevate geometry
      Base::elevate(json);
 
      // Set geometry
      bc_.setGeoMap(Base::geo_);
    }
 
    int num(1), dim(-1), patchIndex(-1);
 
    if (json.contains("data")) {
      if (json["data"].contains("num"))
        num = json["data"]["num"].get<int>();
 
      if (json["data"].contains("dim"))
        dim = json["data"]["dim"].get<int>();
 
      if (json["data"].contains("patch"))
        patchIndex = json["data"]["patch"].get<int>();
    }
 
    // Degree elevate basis of solution space
    if (patchIndex == -1)
      basis_.degreeElevate(num, dim);
    else
      basis_.basis(patchIndex).degreeElevate(num, dim);
 
    // Regenerate solution
    solve();
  }
 
  void increase(const nlohmann::json &json = NULL) override {
 
    bool geometry = true;
 
    if (json.contains("data"))
      if (json["data"].contains("num"))
        geometry = json["data"]["geometry"].get<bool>();
 
    if (geometry) {
      // Increase geometry
      Base::increase(json);
 
      // Set geometry
      bc_.setGeoMap(Base::geo_);
    }
 
    int num(1), dim(-1), patchIndex(-1);
 
    if (json.contains("data")) {
      if (json["data"].contains("num"))
        num = json["data"]["num"].get<int>();
 
      if (json["data"].contains("dim"))
        dim = json["data"]["dim"].get<int>();
 
      if (json["data"].contains("patch"))
        patchIndex = json["data"]["patch"].get<int>();
    }
 
    // Degree increase basis of solution space
    if (patchIndex == -1)
      basis_.degreeIncrease(num, dim);
    else
      basis_.basis(patchIndex).degreeIncrease(num, dim);
 
    // Regenerate solution
    solve();
  }
 
  void refine(const nlohmann::json &json = NULL) override {
 
    bool geometry = true;
 
    if (json.contains("data"))
      if (json["data"].contains("num"))
        geometry = json["data"]["geometry"].get<bool>();
 
    if (geometry) {
      // Refine geometry
      Base::refine(json);
 
      // Set geometry
      bc_.setGeoMap(Base::geo_);
    }
 
    int num(1), dim(-1), patchIndex(-1);
 
    if (json.contains("data")) {
      if (json["data"].contains("num"))
        num = json["data"]["num"].get<int>();
 
      if (json["data"].contains("dim"))
        dim = json["data"]["dim"].get<int>();
 
      if (json["data"].contains("patch"))
        patchIndex = json["data"]["patch"].get<int>();
    }
 
    // Refine basis of solution space
    if (patchIndex == -1)
      basis_.uniformRefine(num, 1, dim);
    else
      basis_.basis(patchIndex).uniformRefine(num, 1, dim);
 
    // Regenerate solution
    solve();
  }
 
  void addPatch(const nlohmann::json &json = NULL) override {
 
    // Add patch from geometry
    Base::addPatch(json);
 
    // Set geometry
    bc_.setGeoMap(Base::geo_);
    
    throw std::runtime_error("Adding patches is not yet implemented in G+Smo");
 
    // Regenerate solution
    solve();
  }
  
  void removePatch(const std::string &patch,
                   const nlohmann::json &json = NULL) override {
 
    // Remove patch from geometry
    Base::removePatch(json);
 
    // Set geometry
    bc_.setGeoMap(Base::geo_);
    
    int patchIndex(-1);
 
    if (json.contains("data")) {      
      if (json["data"].contains("patch"))
        patchIndex = json["data"]["patch"].get<int>();
    }
 
    // if (patchIndex == -1)
    //   throw std::runtime_error("Invalid patch index");
 
    // throw std::runtime_error("Patch removal is not yet implemented in G+Smo");
 
    // Regenerate solution
    solve();
  }
};
 
} // namespace webapp
} // namespace iganet