<commit_before>//
// PROJECT CHRONO - http://projectchrono.org
//
// Copyright (c) 2010-2012 <NAME>
// Copyright (c) 2013 Project Chrono
// 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.
//

#include "chrono/solver/ChSystemDescriptor.h"
#include "chrono/solver/ChConstraintTwoTuplesContactN.h"
#include "chrono/solver/ChConstraintTwoTuplesFrictionT.h"
#include "chrono/core/ChLinkedListMatrix.h"

namespace chrono {

// Register into the object factory, to enable run-time
// dynamic creation and persistence
ChClassRegister<ChSystemDescriptor> a_registration_ChSystemDescriptor;

#define CH_SPINLOCK_HASHSIZE 203

ChSystemDescriptor::ChSystemDescriptor() {
    vconstraints.clear();
    vvariables.clear();
    vstiffness.clear();

    c_a = 1.0;

    n_q = 0;
    n_c = 0;
    freeze_count = false;

    this->num_threads = CHOMPfunctions::GetNumProcs();

    spinlocktable = new ChSpinlock[CH_SPINLOCK_HASHSIZE];
}

ChSystemDescriptor::~ChSystemDescriptor() {
    vconstraints.clear();
    vvariables.clear();
    vstiffness.clear();

    if (spinlocktable)
        delete[] spinlocktable;
    spinlocktable = 0;
}

void ChSystemDescriptor::ComputeFeasabilityViolation(
    double& resulting_maxviolation,  ///< gets the max constraint violation (either bi- and unilateral.)
    double& resulting_feasability    ///< gets the max feasability as max |l*c|, for unilateral only
    ) {
    resulting_maxviolation = 0;
    resulting_feasability = 0;

    for (unsigned int ic = 0; ic < vconstraints.size(); ic++) {
        // the the residual of the constraint..
        double mres_i = vconstraints[ic]->Compute_c_i();

        double candidate_violation = fabs(vconstraints[ic]->Violation(mres_i));

        if (candidate_violation > resulting_maxviolation)
            resulting_maxviolation = candidate_violation;

        if (vconstraints[ic]->IsUnilateral()) {
            double candidate_feas = fabs(mres_i * vconstraints[ic]->Get_l_i());  // =|c*l|
            if (candidate_feas > resulting_feasability)
                resulting_feasability = candidate_feas;
        }
    }
}

int ChSystemDescriptor::CountActiveVariables() {
    if (this->freeze_count)  // optimization, avoid list count all times
        return n_q;

    n_q = 0;
    for (unsigned int iv = 0; iv < vvariables.size(); iv++) {
        if (vvariables[iv]->IsActive()) {
            vvariables[iv]->SetOffset(n_q);  // also store offsets in state and MC matrix
            n_q += vvariables[iv]->Get_ndof();
        }
    }
    return n_q;
}

int ChSystemDescriptor::CountActiveConstraints() {
    if (this->freeze_count)  // optimization, avoid list count all times
        return n_c;

    n_c = 0;
    for (unsigned int ic = 0; ic < vconstraints.size(); ic++) {
        if (vconstraints[ic]->IsActive()) {
            vconstraints[ic]->SetOffset(n_c);  // also store offsets in state and MC matrix
            n_c++;
        }
    }
    return n_c;
}

void ChSystemDescriptor::UpdateCountsAndOffsets() {
    freeze_count = false;
    CountActiveVariables();
    CountActiveConstraints();
    freeze_count = true;
}

void ChSystemDescriptor::ConvertToMatrixForm(ChSparseMatrix* Cq,
                                             ChSparseMatrix* M,
                                             ChSparseMatrix* E,
                                             ChMatrix<>* Fvector,
                                             ChMatrix<>* Bvector,
                                             ChMatrix<>* Frict,
                                             bool only_bilaterals,
                                             bool skip_contacts_uv) {
    std::vector<ChConstraint*>& mconstraints = this->GetConstraintsList();
    std::vector<ChVariables*>& mvariables = this->GetVariablesList();

    // Count bilateral and other constraints.. (if wanted, bilaterals only)

    int mn_c = 0;
    for (unsigned int ic = 0; ic < mconstraints.size(); ic++) {
        if (mconstraints[ic]->IsActive())
            if (!((mconstraints[ic]->GetMode() == CONSTRAINT_FRIC) && only_bilaterals))
                if (!((dynamic_cast<ChConstraintTwoTuplesFrictionTall*>(mconstraints[ic])) && skip_contacts_uv)) {
                    mn_c++;
                }
    }

    // Count active variables, by scanning through all variable blocks,
    // and set offsets

    n_q = this->CountActiveVariables();

    // Reset and resize (if needed) auxiliary vectors

    if (Cq)
        Cq->Reset(mn_c, n_q);
    if (M)
        M->Reset(n_q, n_q);
    if (E)
        E->Reset(mn_c, mn_c);
    if (Fvector)
        Fvector->Reset(n_q, 1);
    if (Bvector)
        Bvector->Reset(mn_c, 1);
    if (Frict)
        Frict->Reset(mn_c, 1);

    // Fills M submasses and 'f' vector,
    // by looping on variables
    int s_q = 0;
    for (unsigned int iv = 0; iv < mvariables.size(); iv++) {
        if (mvariables[iv]->IsActive()) {
            if (M)
                mvariables[iv]->Build_M(*M, s_q, s_q, this->c_a);  // .. fills  M
            if (Fvector)
                Fvector->PasteMatrix(&vvariables[iv]->Get_fb(), s_q, 0);  // .. fills  'f'
            s_q += mvariables[iv]->Get_ndof();
        }
    }

    // If some stiffness / hessian matrix has been added to M ,
    // also add it to the sparse M
    int s_k = 0;
    for (unsigned int ik = 0; ik < this->vstiffness.size(); ik++) {
        this->vstiffness[ik]->Build_K(*M, true);
    }

    // Fills Cq jacobian, E 'compliance' matrix , the 'b' vector and friction coeff.vector,<commit_msg>Optional build of 'mass' matrix H<commit_after>//
// PROJECT CHRONO - http://projectchrono.org
//
// Copyright (c) 2010-2012 <NAME>
// Copyright (c) 2013 Project Chrono
// 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.
//

#include "chrono/solver/ChSystemDescriptor.h"
#include "chrono/solver/ChConstraintTwoTuplesContactN.h"
#include "chrono/solver/ChConstraintTwoTuplesFrictionT.h"
#include "chrono/core/ChLinkedListMatrix.h"

namespace chrono {

// Register into the object factory, to enable run-time
// dynamic creation and persistence
ChClassRegister<ChSystemDescriptor> a_registration_ChSystemDescriptor;

#define CH_SPINLOCK_HASHSIZE 203

ChSystemDescriptor::ChSystemDescriptor() {
    vconstraints.clear();
    vvariables.clear();
    vstiffness.clear();

    c_a = 1.0;

    n_q = 0;
    n_c = 0;
    freeze_count = false;

    this->num_threads = CHOMPfunctions::GetNumProcs();

    spinlocktable = new ChSpinlock[CH_SPINLOCK_HASHSIZE];
}

ChSystemDescriptor::~ChSystemDescriptor() {
    vconstraints.clear();
    vvariables.clear();
    vstiffness.clear();

    if (spinlocktable)
        delete[] spinlocktable;
    spinlocktable = 0;
}

void ChSystemDescriptor::ComputeFeasabilityViolation(
    double& resulting_maxviolation,  ///< gets the max constraint violation (either bi- and unilateral.)
    double& resulting_feasability    ///< gets the max feasability as max |l*c|, for unilateral only
    ) {
    resulting_maxviolation = 0;
    resulting_feasability = 0;

    for (unsigned int ic = 0; ic < vconstraints.size(); ic++) {
        // the the residual of the constraint..
        double mres_i = vconstraints[ic]->Compute_c_i();

        double candidate_violation = fabs(vconstraints[ic]->Violation(mres_i));

        if (candidate_violation > resulting_maxviolation)
            resulting_maxviolation = candidate_violation;

        if (vconstraints[ic]->IsUnilateral()) {
            double candidate_feas = fabs(mres_i * vconstraints[ic]->Get_l_i());  // =|c*l|
            if (candidate_feas > resulting_feasability)
                resulting_feasability = candidate_feas;
        }
    }
}

int ChSystemDescriptor::CountActiveVariables() {
    if (this->freeze_count)  // optimization, avoid list count all times
        return n_q;

    n_q = 0;
    for (unsigned int iv = 0; iv < vvariables.size(); iv++) {
        if (vvariables[iv]->IsActive()) {
            vvariables[iv]->SetOffset(n_q);  // also store offsets in state and MC matrix
            n_q += vvariables[iv]->Get_ndof();
        }
    }
    return n_q;
}

int ChSystemDescriptor::CountActiveConstraints() {
    if (this->freeze_count)  // optimization, avoid list count all times
        return n_c;

    n_c = 0;
    for (unsigned int ic = 0; ic < vconstraints.size(); ic++) {
        if (vconstraints[ic]->IsActive()) {
            vconstraints[ic]->SetOffset(n_c);  // also store offsets in state and MC matrix
            n_c++;
        }
    }
    return n_c;
}

void ChSystemDescriptor::UpdateCountsAndOffsets() {
    freeze_count = false;
    CountActiveVariables();
    CountActiveConstraints();
    freeze_count = true;
}

void ChSystemDescriptor::ConvertToMatrixForm(ChSparseMatrix* Cq,
                                             ChSparseMatrix* M,
                                             ChSparseMatrix* E,
                                             ChMatrix<>* Fvector,
                                             ChMatrix<>* Bvector,
                                             ChMatrix<>* Frict,
                                             bool only_bilaterals,
                                             bool skip_contacts_uv) {
    std::vector<ChConstraint*>& mconstraints = this->GetConstraintsList();
    std::vector<ChVariables*>& mvariables = this->GetVariablesList();

    // Count bilateral and other constraints.. (if wanted, bilaterals only)

    int mn_c = 0;
    for (unsigned int ic = 0; ic < mconstraints.size(); ic++) {
        if (mconstraints[ic]->IsActive())
            if (!((mconstraints[ic]->GetMode() == CONSTRAINT_FRIC) && only_bilaterals))
                if (!((dynamic_cast<ChConstraintTwoTuplesFrictionTall*>(mconstraints[ic])) && skip_contacts_uv)) {
                    mn_c++;
                }
    }

    // Count active variables, by scanning through all variable blocks,
    // and set offsets

    n_q = this->CountActiveVariables();

    // Reset and resize (if needed) auxiliary vectors

    if (Cq)
        Cq->Reset(mn_c, n_q);
    if (M)
        M->Reset(n_q, n_q);
    if (E)
        E->Reset(mn_c, mn_c);
    if (Fvector)
        Fvector->Reset(n_q, 1);
    if (Bvector)
        Bvector->Reset(mn_c, 1);
    if (Frict)
        Frict->Reset(mn_c, 1);

    // Fills M submasses and 'f' vector,
    // by looping on variables
    int s_q = 0;
    for (unsigned int iv = 0; iv < mvariables.size(); iv++) {
        if (mvariables[iv]->IsActive()) {
            if (M)
                mvariables[iv]->Build_M(*M, s_q, s_q, this->c_a);  // .. fills  M
            if (Fvector)
                Fvector->PasteMatrix(&vvariables[iv]->Get_fb(), s_q, 0);  // .. fills  'f'
            s_q += mvariables[iv]->Get_ndof();
        }
    }

    // If some stiffness / hessian matrix has been added to M ,
    // also add it to the sparse M
    int s_k = 0;
    if (H) {
        for (unsigned int ik = 0; ik < this->vstiffness.size(); ik++) {
            this->vstiffness[ik]->Build_K(*H, true);
        }
    }

    // Fills Cq jacobian, E 'compliance' matrix , the 'b' vector and friction coeff.vector,
{'fraction_non_alphanumeric': 0.11981405139900009, 'fraction_numerical': 0.006490658714147882, 'mean_word_length': 2.293471981513576, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 87, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5145693', 'n_tokens_mistral': 3600, 'n_tokens_neox': 3429, 'n_words': 858}