<commit_before>// =============================================================================
// PROJECT CHRONO - http://projectchrono.org
//
// Copyright (c) 2014 projectchrono.org
// All rights reserved.
//
// Use of this source code is governed by a BSD-style license that can be found
// in the LICENSE file at the top level of the distribution and at
// http://projectchrono.org/license-chrono.txt.
//
// =============================================================================
// Authors: <NAME>, <NAME>
// =============================================================================

#include "chrono/fea/ChElementBar.h"

namespace chrono {
namespace fea {

ChElementBar::ChElementBar() : area(0.01 * 0.01), density(1000), E(0.01e9), rdamping(0.01), length(0), mass(0) {
    nodes.resize(2);
}

ChElementBar::~ChElementBar() {}

void ChElementBar::SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA, std::shared_ptr<ChNodeFEAxyz> nodeB) {
    nodes[0] = nodeA;
    nodes[1] = nodeB;
    std::vector<ChVariables*> mvars;
    mvars.push_back(&nodes[0]->Variables());
    mvars.push_back(&nodes[1]->Variables());
    Kmatr.SetVariables(mvars);
}

void ChElementBar::GetStateBlock(ChVectorDynamic<>& mD) {
    mD.setZero(this->GetNdofs());
    mD.segment(0, 3) = this->nodes[0]->GetPos().eigen();
    mD.segment(3, 3) = this->nodes[1]->GetPos().eigen();
}

void ChElementBar::ComputeKRMmatricesGlobal(ChMatrixRef H, double Kfactor, double Rfactor, double Mfactor) {
    assert((H.rows() == 6) && (H.cols() == 6));

    // For K stiffness matrix and R damping matrix:
    // compute stiffness matrix (this is already the explicit
    // formulation of the corotational stiffness matrix in 3D)

    ChVector<> dir = (nodes[1]->GetPos() - nodes[0]->GetPos()).GetNormalized();
    ChVectorN<double, 3> dircolumn = dir.eigen();

    double Kstiffness = ((this->area * this->E) / this->length);
    double Rdamping = this->rdamping * Kstiffness;
    // note that stiffness and damping matrices are the same, so join stuff here
    double commonfactor = Kstiffness * Kfactor + Rdamping * Rfactor;
    ChMatrix33<> submatr = commonfactor * dircolumn * dircolumn.transpose();

    H.block(0,0,3,3) = submatr;
    H.block(3,3,3,3) = submatr;
    H.block(0,3,3,3) = -submatr;
    H.block(3,0,3,3) = -submatr;

    // For M mass matrix, do mass lumping:
    H(0, 0) += Mfactor * mass * 0.5;  // node A x,y,z
    H(1, 1) += Mfactor * mass * 0.5;
    H(2, 2) += Mfactor * mass * 0.5;
    H(3, 3) += Mfactor * mass * 0.5;  // node B x,y,z<commit_msg>Add tangent geometric stiffness for the ChElementBar. Previously only the material tangent stiffness was used, but geometric stiffness can help convergence a lot when studying pre-stressed structures. Also, without it, the eigenvalues (natural frequencies) would be wrong in prestressed structures.<commit_after>// =============================================================================
// PROJECT CHRONO - http://projectchrono.org
//
// Copyright (c) 2014 projectchrono.org
// All rights reserved.
//
// Use of this source code is governed by a BSD-style license that can be found
// in the LICENSE file at the top level of the distribution and at
// http://projectchrono.org/license-chrono.txt.
//
// =============================================================================
// Authors: <NAME>, <NAME>
// =============================================================================

#include "chrono/fea/ChElementBar.h"

namespace chrono {
namespace fea {

ChElementBar::ChElementBar() : area(0.01 * 0.01), density(1000), E(0.01e9), rdamping(0.01), length(0), mass(0) {
    nodes.resize(2);
}

ChElementBar::~ChElementBar() {}

void ChElementBar::SetNodes(std::shared_ptr<ChNodeFEAxyz> nodeA, std::shared_ptr<ChNodeFEAxyz> nodeB) {
    nodes[0] = nodeA;
    nodes[1] = nodeB;
    std::vector<ChVariables*> mvars;
    mvars.push_back(&nodes[0]->Variables());
    mvars.push_back(&nodes[1]->Variables());
    Kmatr.SetVariables(mvars);
}

void ChElementBar::GetStateBlock(ChVectorDynamic<>& mD) {
    mD.setZero(this->GetNdofs());
    mD.segment(0, 3) = this->nodes[0]->GetPos().eigen();
    mD.segment(3, 3) = this->nodes[1]->GetPos().eigen();
}

void ChElementBar::ComputeKRMmatricesGlobal(ChMatrixRef H, double Kfactor, double Rfactor, double Mfactor) {
    assert((H.rows() == 6) && (H.cols() == 6));

    // For K stiffness matrix and R damping matrix:
    // compute stiffness matrix (this is already the explicit
    // formulation of the corotational stiffness matrix in 3D)

    ChVector<> dir = (nodes[1]->GetPos() - nodes[0]->GetPos()).GetNormalized();
    ChVectorN<double, 3> dircolumn = dir.eigen();

    double Kstiffness = ((this->area * this->E) / this->length);
    double Rdamping = this->rdamping * Kstiffness;
    // note that stiffness and damping matrices are the same, so join stuff here
    double commonfactor = Kstiffness * Kfactor + Rdamping * Rfactor;
	ChMatrix33<> V = dircolumn * dircolumn.transpose();
    ChMatrix33<> keV = commonfactor * V;

    H.block(0,0,3,3) = keV;
    H.block(3,3,3,3) = keV;
    H.block(0,3,3,3) = -keV;
    H.block(3,0,3,3) = -keV;

	// add geometric stiffness - in future it might become an option to switch off if not needed.
	// See for ex. http://shodhbhagirathi.iitr.ac.in:8081/jspui/handle/123456789/8433 pag. 14-15
	if (true) {
		double L_ref = (nodes[1]->GetX0() - nodes[0]->GetX0()).Length();
		double L = (nodes[1]->GetPos() - nodes[0]->GetPos()).Length();
		double Kstiffness = ((this->area * this->E) / this->length);
		double internal_Kforce_local = Kstiffness * (L - L_ref);

		ChMatrix33<> kgV = Kfactor * (internal_Kforce_local / L_ref) * (ChMatrix33<>(1) - V);
		H.block(0, 0, 3, 3) += kgV;
		H.block(3, 3, 3, 3) += kgV;
		H.block(0, 3, 3, 3) += -kgV;
		H.block(3, 0, 3, 3) += -kgV;
	}

    // For M mass matrix, do mass lumping:
    H(0, 0) += Mfactor * mass * 0.5;  // node A x,y,z
    H(1, 1) += Mfactor * mass * 0.5;
    H(2, 2) += Mfactor * mass * 0.5;
    H(3, 3) += Mfactor * mass * 0.5;  // node B x,y,z
{'fraction_non_alphanumeric': 0.2310750372701673, 'fraction_numerical': 0.032797747225443104, 'mean_word_length': 3.9451269451269453, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2207691', 'n_tokens_mistral': 2171, 'n_tokens_neox': 2010, 'n_words': 514}