<commit_before>
#include <algorithm>
#include <math.h>
#include "chrono_parallel/physics/ChSystemParallel.h"
#include <chrono_parallel/physics/Ch3DOFContainer.h>
#include "chrono_parallel/ChDataManager.h"
#include "chrono_parallel/constraints/ChConstraintFluidFluidUtils.h"
#include "chrono_parallel/constraints/ChConstraintUtils.h"

#include "chrono_parallel/math/other_types.h"  // for uint, int2, vec3
#include "chrono_parallel/math/real.h"         // for real
#include "chrono_parallel/math/real2.h"        // for real2
#include "chrono_parallel/math/real3.h"        // for real3
#include "chrono_parallel/math/real4.h"        // for quaternion, real4
#include "chrono_parallel/math/matrix.h"       // for quaternion, real4
#include "chrono_parallel/collision/ChCollision.h"

namespace chrono {

using namespace collision;
using namespace geometry;

ChFluidContainer::ChFluidContainer(ChSystemParallelDVI* physics_system) {
    data_manager = physics_system->data_manager;
    data_manager->Add3DOFContainer(this);
    body_offset = 0;
    compliance = 0;
    epsilon = 10e-3;
    tau = 4 * .01;

    rho = 1000;
    mass = 0.037037;

    viscosity = 0;

    artificial_pressure = false;
    artificial_pressure_k = 0.1;
    artificial_pressure_dq = .2 * kernel_radius;
    artificial_pressure_n = 4;
    enable_viscosity = false;
    initialize_mass = false;
}
ChFluidContainer::~ChFluidContainer() {}

void ChFluidContainer::AddBodies(const std::vector<real3>& positions, const std::vector<real3>& velocities) {
    custom_vector<real3>& pos_fluid = data_manager->host_data.pos_3dof;
    custom_vector<real3>& vel_fluid = data_manager->host_data.vel_3dof;

    pos_fluid.insert(pos_fluid.end(), positions.begin(), positions.end());
    vel_fluid.insert(vel_fluid.end(), velocities.begin(), velocities.end());
    // In case the number of velocities provided were not enough, resize to the number of fluid bodies
    vel_fluid.resize(pos_fluid.size());
    data_manager->num_fluid_bodies = pos_fluid.size();
}
void ChFluidContainer::Update(double ChTime) {
    uint num_fluid_bodies = data_manager->num_fluid_bodies;
    uint num_rigid_bodies = data_manager->num_rigid_bodies;
    uint num_shafts = data_manager->num_shafts;
    custom_vector<real3>& pos_fluid = data_manager->host_data.pos_3dof;
    custom_vector<real3>& vel_fluid = data_manager->host_data.vel_3dof;
    real3 g_acc = data_manager->settings.gravity;
    real3 h_gravity = data_manager->settings.step_size * mass * g_acc;
    for (int i = 0; i < num_fluid_bodies; i++) {
        // This was moved to after fluid collision detection
        // real3 vel = vel_fluid[i];
        // data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 0] = vel.x;
        // data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 1] = vel.y;
        // data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 2] = vel.z;

        data_manager->host_data.hf[num_rigid_bodies * 6 + num_shafts + i * 3 + 0] = h_gravity.x;
        data_manager->host_data.hf[num_rigid_bodies * 6 + num_shafts + i * 3 + 1] = h_gravity.y;
        data_manager->host_data.hf[num_rigid_bodies * 6 + num_shafts + i * 3 + 2] = h_gravity.z;
    }
}

void ChFluidContainer::UpdatePosition(double ChTime) {
    uint num_fluid_bodies = data_manager->num_fluid_bodies;
    uint num_rigid_bodies = data_manager->num_rigid_bodies;
    uint num_shafts = data_manager->num_shafts;

    custom_vector<real3>& pos_fluid = data_manager->host_data.pos_3dof;
    custom_vector<real3>& vel_fluid = data_manager->host_data.vel_3dof;
    custom_vector<real3>& sorted_pos_fluid = data_manager->host_data.sorted_pos_3dof;
#pragma omp parallel for
    for (int i = 0; i < num_fluid_bodies; i++) {
        real3 vel;
        int original_index = data_manager->host_data.particle_indices_3dof[i];
        // these are sorted so we have to unsort them
        vel.x = data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 0];
        vel.y = data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 1];
        vel.z = data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 2];

        real speed = Length(vel);
        if (speed > max_velocity) {
            vel = vel * max_velocity / speed;
        }
        vel_fluid[original_index] = vel;
        pos_fluid[original_index] += vel * data_manager->settings.step_size;
        sorted_pos_fluid[i] = pos_fluid[original_index];
    }
    new_pos = sorted_pos_fluid;
    if (data_manager->num_fluid_bodies != 0) {
        data_manager->narrowphase->DispatchRigidFluid();
    }

    custom_vector<real3>& cpta = data_manager->host_data.cpta_rigid_fluid;
    custom_vector<real3>& norm = data_manager->host_data.norm_rigid_fluid;
    custom_vector<real>& dpth = data_manager->host_data.dpth_rigid_fluid;
    custom_vector<int>& neighbor_rigid_fluid = data_manager->host_data.neighbor_rigid_fluid;
    custom_vector<int>& contact_counts = data_manager->host_data.c_counts_rigid_fluid;
    // This treats all rigid neighbors as fixed. This correction should usually be pretty small if the timestep isnt
    // too large.

    if (data_manager->num_rigid_fluid_contacts > 0) {
#pragma omp parallel for
        for (int p = 0; p < num_fluid_bodies; p++) {
            int start = contact_counts[p];
            int end = contact_counts[p + 1];
            real3 delta = real3(0);
            real weight = 0;
            for (int index = start; index < end; index++) {
                int i = index - start;
                int rigid = neighbor_rigid_fluid[p * max_rigid_neighbors + i];
                // if (data_manager->host_data.active_rigid[rigid] == false) {
                real3 U = norm[p * max_rigid_neighbors + i];
                real depth = dpth[p * max_rigid_neighbors + i];
                if (depth < 0) {
                    real w = 1.0;  // mass / (mass + data_manager->host_data.mass_rigid[rigid]);
                    delta -= w * depth * U;
                    weight++;
                }
                //}
            }
            if (weight > 0) {
                new_pos[p] = new_pos[p] + delta / weight;
            }
        }

#pragma omp parallel for
        for (int p = 0; p < num_fluid_bodies; p++) {
            // real3 vnew = (new_pos[p] - pos_fluid[p]) / data_manager->settings.step_size;
            int original_index = data_manager->host_data.particle_indices_3dof[p];
            pos_fluid[original_index] = new_pos[p];
        }
    }
}

int ChFluidContainer::GetNumConstraints() {
    int num_fluid_fluid = data_manager->num_fluid_bodies + data_manager->num_rigid_fluid_contacts * 3;<commit_msg>Disable position based correction for fluid<commit_after>
#include <algorithm>
#include <math.h>
#include "chrono_parallel/physics/ChSystemParallel.h"
#include <chrono_parallel/physics/Ch3DOFContainer.h>
#include "chrono_parallel/ChDataManager.h"
#include "chrono_parallel/constraints/ChConstraintFluidFluidUtils.h"
#include "chrono_parallel/constraints/ChConstraintUtils.h"

#include "chrono_parallel/math/other_types.h"  // for uint, int2, vec3
#include "chrono_parallel/math/real.h"         // for real
#include "chrono_parallel/math/real2.h"        // for real2
#include "chrono_parallel/math/real3.h"        // for real3
#include "chrono_parallel/math/real4.h"        // for quaternion, real4
#include "chrono_parallel/math/matrix.h"       // for quaternion, real4
#include "chrono_parallel/collision/ChCollision.h"

namespace chrono {

using namespace collision;
using namespace geometry;

ChFluidContainer::ChFluidContainer(ChSystemParallelDVI* physics_system) {
    data_manager = physics_system->data_manager;
    data_manager->Add3DOFContainer(this);
    body_offset = 0;
    compliance = 0;
    epsilon = 10e-3;
    tau = 4 * .01;

    rho = 1000;
    mass = 0.037037;

    viscosity = 0;

    artificial_pressure = false;
    artificial_pressure_k = 0.1;
    artificial_pressure_dq = .2 * kernel_radius;
    artificial_pressure_n = 4;
    enable_viscosity = false;
    initialize_mass = false;
}
ChFluidContainer::~ChFluidContainer() {}

void ChFluidContainer::AddBodies(const std::vector<real3>& positions, const std::vector<real3>& velocities) {
    custom_vector<real3>& pos_fluid = data_manager->host_data.pos_3dof;
    custom_vector<real3>& vel_fluid = data_manager->host_data.vel_3dof;

    pos_fluid.insert(pos_fluid.end(), positions.begin(), positions.end());
    vel_fluid.insert(vel_fluid.end(), velocities.begin(), velocities.end());
    // In case the number of velocities provided were not enough, resize to the number of fluid bodies
    vel_fluid.resize(pos_fluid.size());
    data_manager->num_fluid_bodies = pos_fluid.size();
}
void ChFluidContainer::Update(double ChTime) {
    uint num_fluid_bodies = data_manager->num_fluid_bodies;
    uint num_rigid_bodies = data_manager->num_rigid_bodies;
    uint num_shafts = data_manager->num_shafts;
    custom_vector<real3>& pos_fluid = data_manager->host_data.pos_3dof;
    custom_vector<real3>& vel_fluid = data_manager->host_data.vel_3dof;
    real3 g_acc = data_manager->settings.gravity;
    real3 h_gravity = data_manager->settings.step_size * mass * g_acc;
    for (int i = 0; i < num_fluid_bodies; i++) {
        // This was moved to after fluid collision detection
        // real3 vel = vel_fluid[i];
        // data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 0] = vel.x;
        // data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 1] = vel.y;
        // data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 2] = vel.z;

        data_manager->host_data.hf[num_rigid_bodies * 6 + num_shafts + i * 3 + 0] = h_gravity.x;
        data_manager->host_data.hf[num_rigid_bodies * 6 + num_shafts + i * 3 + 1] = h_gravity.y;
        data_manager->host_data.hf[num_rigid_bodies * 6 + num_shafts + i * 3 + 2] = h_gravity.z;
    }
}

void ChFluidContainer::UpdatePosition(double ChTime) {
    uint num_fluid_bodies = data_manager->num_fluid_bodies;
    uint num_rigid_bodies = data_manager->num_rigid_bodies;
    uint num_shafts = data_manager->num_shafts;

    custom_vector<real3>& pos_fluid = data_manager->host_data.pos_3dof;
    custom_vector<real3>& vel_fluid = data_manager->host_data.vel_3dof;
    custom_vector<real3>& sorted_pos_fluid = data_manager->host_data.sorted_pos_3dof;
#pragma omp parallel for
    for (int i = 0; i < num_fluid_bodies; i++) {
        real3 vel;
        int original_index = data_manager->host_data.particle_indices_3dof[i];
        // these are sorted so we have to unsort them
        vel.x = data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 0];
        vel.y = data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 1];
        vel.z = data_manager->host_data.v[num_rigid_bodies * 6 + num_shafts + i * 3 + 2];

        real speed = Length(vel);
        if (speed > max_velocity) {
            vel = vel * max_velocity / speed;
        }
        vel_fluid[original_index] = vel;
        pos_fluid[original_index] += vel * data_manager->settings.step_size;
        // sorted_pos_fluid[i] = pos_fluid[original_index];
    }
    //    new_pos = sorted_pos_fluid;
    //    if (data_manager->num_fluid_bodies != 0) {
    //        data_manager->narrowphase->DispatchRigidFluid();
    //    }
    //
    //    custom_vector<real3>& cpta = data_manager->host_data.cpta_rigid_fluid;
    //    custom_vector<real3>& norm = data_manager->host_data.norm_rigid_fluid;
    //    custom_vector<real>& dpth = data_manager->host_data.dpth_rigid_fluid;
    //    custom_vector<int>& neighbor_rigid_fluid = data_manager->host_data.neighbor_rigid_fluid;
    //    custom_vector<int>& contact_counts = data_manager->host_data.c_counts_rigid_fluid;
    //    // This treats all rigid neighbors as fixed. This correction should usually be pretty small if the timestep
    //    isnt
    //    // too large.
    //
    //    if (data_manager->num_rigid_fluid_contacts > 0) {
    //#pragma omp parallel for
    //        for (int p = 0; p < num_fluid_bodies; p++) {
    //            int start = contact_counts[p];
    //            int end = contact_counts[p + 1];
    //            real3 delta = real3(0);
    //            real weight = 0;
    //            for (int index = start; index < end; index++) {
    //                int i = index - start;
    //                int rigid = neighbor_rigid_fluid[p * max_rigid_neighbors + i];
    //                // if (data_manager->host_data.active_rigid[rigid] == false) {
    //                real3 U = norm[p * max_rigid_neighbors + i];
    //                real depth = dpth[p * max_rigid_neighbors + i];
    //                if (depth < 0) {
    //                    real w = 1.0;  // mass / (mass + data_manager->host_data.mass_rigid[rigid]);
    //                    delta -= w * depth * U;
    //                    weight++;
    //                }
    //                //}
    //            }
    //            if (weight > 0) {
    //                new_pos[p] = new_pos[p] + delta / weight;
    //            }
    //        }
    //
    //#pragma omp parallel for
    //        for (int p = 0; p < num_fluid_bodies; p++) {
    //            // real3 vnew = (new_pos[p] - pos_fluid[p]) / data_manager->settings.step_size;
    //            int original_index = data_manager->host_data.particle_indices_3dof[p];
    //            pos_fluid[original_index] = new_pos[p];
    //        }
    //    }
}

int ChFluidContainer::GetNumConstraints() {
    int num_fluid_fluid = data_manager->num_fluid_bodies + data_manager->num_rigid_fluid_contacts * 3;
{'fraction_non_alphanumeric': 0.1170018281535649, 'fraction_numerical': 0.015210237659963438, 'mean_word_length': 2.8383384788099915, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 131, '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': '5005043', 'n_tokens_mistral': 4971, 'n_tokens_neox': 4634, 'n_words': 888}