Commit 56cca5e4 authored by Guido Giuntoli's avatar Guido Giuntoli

Merge branch 'mix-laws'

parents 464ecc0a 03e8e188
......@@ -30,6 +30,7 @@
#include <iomanip>
#include <string>
#include <sstream>
#include <map>
#include <cmath>
#include <cassert>
......@@ -131,19 +132,39 @@ typedef struct {
enum {
MIC_HOMOGENEOUS,
MIC_SPHERE,
MIC_LAYER_Y,
MIC_CILI_FIB_Z,
MIC_CILI_FIB_XZ,
MIC_QUAD_FIB_XYZ,
MIC_QUAD_FIB_XZ,
MIC_SPHERE,
MIC_LAYER_Y,
MIC_CILI_FIB_X,
MIC_CILI_FIB_Z,
MIC_CILI_FIB_XZ,
MIC_QUAD_FIB_XYZ,
MIC_QUAD_FIB_XZ,
MIC_QUAD_FIB_XZ_BROKEN_X,
MIC3D_SPHERES,
MIC3D_8,
MIC3D_FIBS_20_ORDER,
MIC3D_FIBS_20_DISORDER
};
static map<int, std::string> micro_names = {
{MIC_HOMOGENEOUS, "MIC_HOMOGENEOUS"},
{MIC_SPHERE, "MIC_SPHERE"},
{MIC_LAYER_Y, "MIC_LAYER_Y"},
{MIC_CILI_FIB_X, "MIC_CILI_FIB_X"},
{MIC_CILI_FIB_Z, "MIC_CILI_FIB_Z"},
{MIC_CILI_FIB_XZ, "MIC_CILI_FIB_XZ"},
{MIC_QUAD_FIB_XYZ, "MIC_QUAD_FIB_XYZ"},
{MIC_QUAD_FIB_XZ, "MIC_QUAD_FIB_XZ"},
{MIC_QUAD_FIB_XZ_BROKEN_X, "MIC_QUAD_FIB_XZ_BROKEN_X"},
{MIC3D_SPHERES, "MIC3D_SPHERES"},
{MIC3D_8, "MIC3D_8"},
{MIC3D_FIBS_20_ORDER, "MIC3D_FIBS_20_ORDER"},
{MIC3D_FIBS_20_DISORDER, "MIC3D_FIBS_20_DISORDER"}
};
/*
* MIC_SPHERES : (2 materials) One sphere in the middle
*
* MIC_HOMOGENEOUS : Only one material (mat[0])
*
* MIC3D_SPHERES : (2 materials) Random spheres.
......@@ -151,15 +172,26 @@ enum {
* MIC3D_8 : (3 materiales) 2 cilinders at 90 deg with a layer around the
* perimeter and a flat layer between the fibers.
*
* MIC3D_FIBS_20_ORDER: (2 materiales) 20 fibers in oriented in X direction.
*
* MIC3D_FIBS_20_DISORDER: (2 materiales) 20 fibers in random directions.
*
*/
enum {
LINEAR,
ONE_WAY,
FULL,
RULE_MIXTURE_1
FE_LINEAR,
FE_ONE_WAY,
FE_FULL,
MIX_RULE_CHAMIS
};
static map<int, int> gp_counter = {
{FE_LINEAR, 0},
{FE_ONE_WAY, 0},
{FE_FULL, 0},
{MIX_RULE_CHAMIS, 0}
};
......@@ -193,7 +225,7 @@ class micropp {
double geo_params[num_geo_params];
material_t *material_list[MAX_MATERIALS];
double ctan_lin[nvoi * nvoi];
double ctan_lin_fe[nvoi * nvoi];
int *elem_type;
double *elem_stress;
......@@ -216,11 +248,6 @@ class micropp {
const bool lin_stress;
/* Number of micro-problems depending on the type */
int num_no_coupling = 0;
int num_one_way = 0;
int num_full = 0;
/* Linear jacobian for optimization */
bool use_A0;
int its_with_A0;
......@@ -235,17 +262,25 @@ class micropp {
int log_id = 0;
ofstream ofstream_log;
/* GPU number of device selection */
/* GPU number for device selection */
int gpu_id = 0;
/* Private function members */
/*
* Linear homogenizations
* Applies to FE RVE model and Mixture rules
*
*/
void homogenize_linear(gp_t<tdim> *gp_ptr);
void homogenize_non_linear(gp_t<tdim> *gp_ptr);
void homogenize_rule_mixture_1(gp_t<tdim> *gp_ptr);
void calc_ctan_lin();
/* FE-based homogenizations */
void homogenize_fe_one_way(gp_t<tdim> *gp_ptr);
void homogenize_fe_full(gp_t<tdim> *gp_ptr);
void calc_ctan_lin_fe_models();
void calc_ctan_lin_mix_rule_Chamis(double ctan[nvoi * nvoi]);
material_t *get_material(const int e) const;
......@@ -333,6 +368,8 @@ class micropp {
void homogenize();
void homogenize_linear();
/* Extras */
int is_non_linear(const int gp_id) const;
......
......@@ -51,6 +51,8 @@ extern "C" {
void micropp3_homogenize(struct micropp3 *self);
void micropp3_homogenize_linear(struct micropp3 *self);
void micropp3_update_vars(struct micropp3 *self);
bool micropp3_is_non_linear(const struct micropp3 *self,
......
......@@ -114,7 +114,7 @@ inline bool point_inside_sphere(const double center[3], const double radius,
const double point[3])
{
/*
* Returns <true> if <point> is inside the sphere with <center> and
* Returns <true> if <point> is inside the sphere with <center> and
* <radius>. Returns <false> otherwise.
*/
......@@ -131,7 +131,7 @@ inline bool point_inside_cilinder_inf(const double dir[3], const double center[3
const double radius, const double point[3])
{
/*
* Returns <true> if <point> is inside the infinite cilinder with
* Returns <true> if <point> is inside the infinite cilinder with
* <direction>, <center> and <radius>. Returns <false> otherwise.
*/
......@@ -151,4 +151,27 @@ inline bool point_inside_cilinder_inf(const double dir[3], const double center[3
}
inline double invert_3x3(const double mat[3][3], double mat_inv[3][3])
{
const double det =
mat[0][0] * (mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2]) -
mat[0][1] * (mat[1][0] * mat[2][2] - mat[2][0] * mat[1][2]) +
mat[0][2] * (mat[1][0] * mat[2][1] - mat[2][0] * mat[1][1]);
mat_inv[0][0] = +(mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2]) / det;
mat_inv[0][1] = -(mat[1][0] * mat[2][2] - mat[2][0] * mat[1][2]) / det;
mat_inv[0][2] = +(mat[1][0] * mat[2][1] - mat[2][0] * mat[1][1]) / det;
mat_inv[1][0] = -(mat[0][1] * mat[2][2] - mat[2][1] * mat[0][2]) / det;
mat_inv[1][1] = +(mat[0][0] * mat[2][2] - mat[2][0] * mat[0][2]) / det;
mat_inv[1][2] = -(mat[0][0] * mat[2][1] - mat[2][0] * mat[0][1]) / det;
mat_inv[2][0] = +(mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2]) / det;
mat_inv[2][1] = -(mat[0][0] * mat[1][2] - mat[1][0] * mat[0][2]) / det;
mat_inv[2][2] = +(mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]) / det;
return det;
}
#endif
......@@ -21,9 +21,11 @@
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*/
#include <cmath>
#include <cassert>
#include "instrument.hpp"
#include "micro.hpp"
......@@ -56,7 +58,7 @@ void micropp<tdim>::get_ctan(const int gp_id, double *ctan) const
template <int tdim>
void micropp<tdim>::homogenize()
void micropp<tdim>::homogenize_linear()
{
INST_START;
......@@ -65,22 +67,47 @@ void micropp<tdim>::homogenize()
gp_t<tdim> *gp_ptr = &gp_list[igp];
if (gp_ptr->coupling == RULE_MIXTURE_1) {
/*
* Computational cheap calculation
* stress = ctan_lin * strain
*/
homogenize_linear(gp_ptr);
}
}
template <int tdim>
void micropp<tdim>::homogenize()
{
INST_START;
#pragma omp parallel for schedule(dynamic,1)
for (int igp = 0; igp < ngp; ++igp) {
homogenize_rule_mixture_1(gp_ptr);
gp_t<tdim> *gp_ptr = &gp_list[igp];
} else if (gp_ptr->coupling == LINEAR) {
if (gp_ptr->coupling == FE_LINEAR ||
gp_ptr->coupling == MIX_RULE_CHAMIS) {
/*
* Computational cheap calculation
* stress = ctan_lin * strain
*
* All mixture rules are linear in Micropp
* so the homogenization of the stress tensor
* is this simple and cheap procedure.
*/
homogenize_linear(gp_ptr);
} else {
} else if (gp_ptr->coupling == FE_ONE_WAY) {
homogenize_fe_one_way(gp_ptr);
homogenize_non_linear(gp_ptr);
} else if (gp_ptr->coupling == FE_FULL) {
homogenize_fe_full(gp_ptr);
}
}
......@@ -97,15 +124,98 @@ void micropp<tdim>::homogenize_linear(gp_t<tdim> * gp_ptr)
memset (gp_ptr->stress, 0.0, nvoi * sizeof(double));
for (int i = 0; i < nvoi; ++i) {
for (int j = 0; j < nvoi; ++j) {
gp_ptr->stress[i] +=
ctan_lin[i * nvoi + j] * gp_ptr->strain[j];
gp_ptr->stress[i] += gp_ptr->ctan[i * nvoi + j] * gp_ptr->strain[j];
}
}
}
template<int tdim>
void micropp<tdim>::homogenize_fe_one_way(gp_t<tdim> * gp_ptr)
{
ell_matrix A; // Jacobian
const int ns[3] = { nx, ny, nz };
ell_init(&A, dim, dim, ns, CG_ABS_TOL, CG_REL_TOL, CG_MAX_ITS);
double *b = (double *) calloc(nndim, sizeof(double));
double *du = (double *) calloc(nndim, sizeof(double));
double *u = (double *) calloc(nndim, sizeof(double));
double *vars_new_aux = (double *) calloc(nvars, sizeof(double));
double *vars_new = (gp_ptr->allocated) ? gp_ptr->vars_k : vars_new_aux;
gp_ptr->cost = 0;
gp_ptr->subiterated = false;
// SIGMA 1 Newton-Raphson
memcpy(u, gp_ptr->u_n, nndim * sizeof(double));
newton_t newton = newton_raphson(&A, b, u, du, gp_ptr->strain, gp_ptr->vars_n);
memcpy(gp_ptr->u_k, u, nndim * sizeof(double));
gp_ptr->cost += newton.solver_its;
gp_ptr->converged = newton.converged;
/*
* In case it has not converged do the sub-iterations
*/
if (gp_ptr->converged == false && subiterations == true) {
double eps_sub[nvoi], deps_sub[nvoi];
memcpy(u, gp_ptr->u_n, nndim * sizeof(double));
memcpy(eps_sub, gp_ptr->strain_old, nvoi * sizeof(double));
gp_ptr->subiterated = true;
for (int i = 0; i < nvoi; ++i)
deps_sub[i] = (gp_ptr->strain[i] - gp_ptr->strain_old[i]) / nsubiterations;
for (int its = 0; its < nsubiterations; ++its) {
for (int j = 0; j < nvoi; ++j) {
eps_sub[j] += deps_sub[j];
}
newton = newton_raphson(&A, b, u, du, eps_sub, gp_ptr->vars_n);
gp_ptr->cost += newton.solver_its;
}
gp_ptr->converged = newton.converged;
memcpy(gp_ptr->u_k, u, nndim * sizeof(double));
}
if (lin_stress) {
memset (gp_ptr->stress, 0.0, nvoi * sizeof(double));
for (int i = 0; i < nvoi; ++i) {
for (int j = 0; j < nvoi; ++j) {
gp_ptr->stress[i] += gp_ptr->ctan[i * nvoi + j] * gp_ptr->strain[j];
}
}
} else {
calc_ave_stress(gp_ptr->u_k, gp_ptr->stress, gp_ptr->vars_n);
}
// Updates <vars_new>
bool non_linear = calc_vars_new(gp_ptr->u_k, gp_ptr->vars_n, vars_new);
if (non_linear == true) {
if (gp_ptr->allocated == false) {
gp_ptr->allocate();
memcpy(gp_ptr->vars_k, vars_new, nvars * sizeof(double));
}
}
ell_free(&A);
free(b);
free(u);
free(du);
free(vars_new_aux);
}
template<int tdim>
void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
void micropp<tdim>::homogenize_fe_full(gp_t<tdim> * gp_ptr)
{
ell_matrix A; // Jacobian
......@@ -124,8 +234,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
// SIGMA 1 Newton-Raphson
memcpy(u, gp_ptr->u_n, nndim * sizeof(double));
newton_t newton = newton_raphson(&A, b, u, du, gp_ptr->strain,
gp_ptr->vars_n);
newton_t newton = newton_raphson(&A, b, u, du, gp_ptr->strain, gp_ptr->vars_n);
memcpy(gp_ptr->u_k, u, nndim * sizeof(double));
gp_ptr->cost += newton.solver_its;
......@@ -162,8 +271,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
memset (gp_ptr->stress, 0.0, nvoi * sizeof(double));
for (int i = 0; i < nvoi; ++i) {
for (int j = 0; j < nvoi; ++j) {
gp_ptr->stress[i] += ctan_lin[i * nvoi + j] *
gp_ptr->strain[j];
gp_ptr->stress[i] += gp_ptr->ctan[i * nvoi + j] * gp_ptr->strain[j];
}
}
......@@ -181,7 +289,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
}
}
if (gp_ptr->allocated && gp_ptr->coupling == FULL) {
if (gp_ptr->allocated) {
// CTAN 3/6 Newton-Raphsons in 2D/3D
double eps_1[6], sig_0[6], sig_1[6];
......@@ -201,8 +309,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
calc_ave_stress(u, sig_1, gp_ptr->vars_n);
for (int v = 0; v < nvoi; ++v)
gp_ptr->ctan[v * nvoi + i] =
(sig_1[v] - sig_0[v]) / D_EPS_CTAN_AVE;
gp_ptr->ctan[v * nvoi + i] = (sig_1[v] - sig_0[v]) / D_EPS_CTAN_AVE;
}
}
......@@ -215,13 +322,6 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
}
template<int tdim>
void micropp<tdim>::homogenize_rule_mixture_1(gp_t<tdim> * gp_ptr)
{
}
template <int tdim>
void micropp<tdim>::update_vars()
{
......
......@@ -91,6 +91,12 @@ module libmicropp
type(micropp3), intent(inout) :: this
end subroutine micropp3_homogenize
subroutine micropp3_homogenize_linear(this) bind(C)
import micropp3
implicit none
type(micropp3), intent(inout) :: this
end subroutine micropp3_homogenize_linear
logical(c_bool) function micropp3_is_non_linear(this, gp_id) bind(C)
use, intrinsic :: iso_c_binding, only: c_bool, c_int
import micropp3
......
This diff is collapsed.
......@@ -53,6 +53,8 @@ extern "C" {
params.write_log = true;
self->ptr = new micropp<3>(params);
delete [] params.coupling;
}
void micropp3_free(micropp3 *self)
......@@ -88,6 +90,12 @@ extern "C" {
ptr->homogenize();
}
void micropp3_homogenize_linear(micropp3 *self)
{
micropp<3> *ptr = (micropp<3> *) self->ptr;
ptr->homogenize_linear();
}
int micropp3_get_cost(const micropp3 *self, int gp_id)
{
micropp<3> *ptr = (micropp<3> *) self->ptr;
......
......@@ -25,6 +25,7 @@ set(testsources
benchmark-mic-1.cpp
benchmark-mic-2.cpp
benchmark-mic-3.cpp
benchmark-mic-4.cpp
)
# Iterate over the list above
......
/*
* This is a test example for MicroPP: a finite element library
* to solve microstructural problems for composite materials.
*
* Copyright (C) - 2018 - Jimmy Aguilar Mena <kratsbinovish@gmail.com>
* Guido Giuntoli <gagiuntoli@gmail.com>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <iomanip>
#include <cstring>
#include <cassert>
#include <chrono>
#include <bits/stdc++.h>
#include "micro.hpp"
using namespace std;
using namespace std::chrono;
double eps_vs_t(double time, double t_final) {
const double eps_max = 1.0e-1;
double eps = eps_max * time;
return eps;
}
int main(int argc, char **argv)
{
if (argc < 2) {
cerr << "Usage: " << argv[0]
<< " n [print=0(def)|1] [steps] [fe_lin=0(def)|fe_one_way|fe_full|rule_mix]" << endl;
return 1;
}
const int n = atoi(argv[1]);
const int print = (argc > 2 ? atoi(argv[2]) : 0); // Optional value
const int time_steps = (argc > 3 ? atoi(argv[3]) : 10); // Optional value
const int coupling = (argc > 4 ? atoi(argv[4]) : 0); // Optional value
const int dir = 0;
const double t_final = 0.15;
const double dt = t_final / time_steps;
double time = 0.0;
assert(n > 1 && 0 <= print && print < 2 && time_steps >= 0 &&
0 <= coupling && coupling <= MIX_RULE_CHAMIS);
micropp_params_t mic_params;
mic_params.ngp = 1;
mic_params.coupling = new int[1];
mic_params.coupling[0] = coupling;
mic_params.size[0] = n;
mic_params.size[1] = n;
mic_params.size[2] = n;
mic_params.type = MIC3D_FIBS_20_ORDER;
mic_params.geo_params[0] = 0.1;
mic_params.geo_params[1] = 0.02;
mic_params.geo_params[2] = 0.01;
material_set(&mic_params.materials[0], 0, 1.0e7, 0.3, 0.0, 0.0, 0.0);
material_set(&mic_params.materials[1], 0, 3.0e7, 0.3, 0.0, 0.0, 0.0);
material_set(&mic_params.materials[2], 0, 3.0e7, 0.3, 0.0, 0.0, 0.0);
mic_params.lin_stress = true;
//mic_params.print();
micropp<3> micro(mic_params);
micro.print_info();
delete [] mic_params.coupling;
if (print) {
char filename[128];
snprintf(filename, 128, "micropp_%d", 0);
micro.output (0, filename);
}
ofstream file;
file.open("result.dat");
file << scientific << setw(14);
auto start = high_resolution_clock::now();
double sig[6], ctan[36];
double eps[6] = { 0. };
cout << scientific;
for (int t = 0; t < time_steps; ++t) {
cout << "Time step = " << t << endl;
eps[dir] = eps_vs_t(time, t_final);
micro.set_strain(0, eps);
cout << "eps = ";
for (int i = 0; i < 6; ++i) {
cout << eps[i] << " ";
}
cout << endl;
cout << "Homogenizing ..." << endl;
micro.homogenize();
int non_linear = micro.is_non_linear(0);
int cost = micro.get_cost(0);
bool has_converged = micro.has_converged(0);
cout << "NL = " << non_linear << endl
<< "Cost = " << cost << endl
<< "Converged = " << has_converged << endl;
cout << "sig = ";
micro.get_stress(0, sig);
for (int i = 0; i < 6; ++i) {
cout << sig[i] << "\t";
}
cout << endl;
micro.get_ctan(0, ctan);
cout << "ctan = " << endl;
for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) {
cout << ctan[i * 6 + j] << "\t";
}
cout << endl;
}
cout << endl;
micro.update_vars();
file << eps[dir] << "\t"
<< sig[dir] << "\t" << endl;
if (print) {
char filename[128];
snprintf(filename, 128, "micropp_%d", t);
micro.output (0, filename);
}
time += dt;
}
auto stop = high_resolution_clock::now();
auto duration = duration_cast<milliseconds>(stop - start);
cout << "time = " << duration.count() << " ms" << endl;
return 0;
}
......@@ -41,7 +41,7 @@ int main (int argc, char *argv[])
const int n = atoi(argv[1]);
const int mic_selected = atoi(argv[2]);
if (mic_selected < 0 || mic_selected > MIC3D_8) {
if (mic_selected < 0 || mic_selected > MIC3D_FIBS_20_DISORDER) {
cerr << "<mic_type = [0 ... 8]>" << endl;
exit(1);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment