diff --git a/VERSIONS b/VERSIONS index a9ceba65..a52ffa4c 100644 --- a/VERSIONS +++ b/VERSIONS @@ -1,6 +1,6 @@ -New in 4.8.1-2 (24 Jan 2024): - +New in 4.8.1 (02 Feb 2024): - general: fixed compatibility with Gmsh (version 4.12.2). +- module -S: added -resmesh orifield,orifieldn. New in 4.8.0 (09 Jan 2024): - module -T: fixed weibull distribution, made minor fixes. @@ -10,7 +10,7 @@ New in 4.8.0 (09 Jan 2024): - module -V: fixed -datanodecoofact. - module -S: changed simulation.config onto simulation.cfg. - documentation: made minor fix. -- general: fixed compatibility with Gmsh (version 4.12.2). +- general: cleaned tests. * Incompatible changes: In -M, changed default value of -order to 2. diff --git a/doc/conf.py b/doc/conf.py index ab4dba6b..faad9dc0 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -11,8 +11,8 @@ import sphinx_rtd_theme project = u'Neper' -version = u'4.8.1-2' -release = u'4.8.1-2' +version = u'4.8.1' +release = u'4.8.1' author = u'Romain Quey' copyright = u'Romain Quey' language = 'en' diff --git a/doc/exprskeys.rst b/doc/exprskeys.rst index 48aad6fb..9bbf915f 100644 --- a/doc/exprskeys.rst +++ b/doc/exprskeys.rst @@ -582,17 +582,19 @@ Available results / keys for nodes are the following: Available results / keys for elements sets are the following: -========================================== ================================================================ ================================== -**Key** **Descriptor** **Apply to** -:data:`ori` average orientation elset, mesh -:data:`gos` grain orientation spread [#gos]_ elset -:data:`anisogos` grain orientation spread computed from :data:`oridisanisoangles` elset -:data:`oridisanisoangles` orientation distribution principal angles elset, mesh -:data:`oridisanisoaxes` orientation distribution principal axes elset, mesh -:data:`oridisanisofact` orientation distribution factor elset, mesh -:data:`odf(=,...)` ODF defined at elements of orientation space (see also below) tess, tesr, mesh, cell, elt, elset -:data:`odfn(=,...)` ODF defined at nodes of orientation space (see also below) tess, tesr, mesh -========================================== ================================================================ ================================== +========================================== ========================================================================== =================================== +**Key** **Descriptor** **Apply to** +:data:`ori` average orientation elset, mesh +:data:`gos` grain orientation spread [#gos]_ elset +:data:`anisogos` grain orientation spread computed from :data:`oridisanisoangles` elset +:data:`oridisanisoangles` orientation distribution principal angles elset, mesh +:data:`oridisanisoaxes` orientation distribution principal axes elset, mesh +:data:`oridisanisofact` orientation distribution factor elset, mesh +:data:`odf(=,...)` ODF defined at elements of orientation space (see also below) tess, tesr, mesh, cell, elt, elset +:data:`odfn(=,...)` ODF defined at nodes of orientation space (see also below) tess, tesr, mesh +:data:`orifield(var=,...)` :data:`` field defined at elements of orientation space (see below) mesh +:data:`orifieldn(var=,...)` :data:`` field defined at nodes of orientation space (see below) mesh +========================================== ========================================================================== =================================== The ODF (:data:`odf` or :data:`odfn`) of a tessellation or mesh is computed over orientation space (provided using :option:`-orispace`) from the orientations of the (tessellation) cells or (mesh) elsets. The (optional) parameters are: @@ -604,6 +606,15 @@ The ODF (:data:`odf` or :data:`odfn`) of a tessellation or mesh is computed over For a cell, element or elset, :data:`odf` returns the value of the ODF of the tessellation or mesh at the corresponding orientation (and simulation step). +The :data:`orifield` and :data:`orifieldn` of a mesh is computed over orientation space (provided using :option:`-orispace`) from the values of the (mesh) elsets. The mandatory parameter is: + +- :data:`var`: the variable, which must be defined for elsets (i.e., have its files in the simulation directory); + +and the optional parameters are: + +- :data:`theta`: the standard deviation of the kernel (in degrees); +- :data:`weight`: the weight of an elset, which can be a real value or an expression based on the :ref:`mesh_keys` (for elsets) -- by default, the volumes of the elsets are used. + .. _rotations_and_orientations: Rotations and Orientations diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3cd21130..18be53c9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,7 +7,7 @@ if(POLICY CMP0077) cmake_policy(SET CMP0077 NEW) endif() -set(NEPER_VERSION \"4.8.1-2\") +set(NEPER_VERSION \"4.8.1\") project(neper) if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 4.8.1) diff --git a/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin2.c b/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin2.c index 94f28418..fbf88122 100644 --- a/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin2.c +++ b/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin2.c @@ -64,6 +64,14 @@ nes_pproc_entity_builtin_elsets (struct SIM *pSim, struct TESS *pTess, &SimRes); } + else if (!strncmp (res, "orifield", 8)) + { + // if mesh, we compute the odf over orientation space + if (!strcmp (entity, "mesh")) + nes_pproc_entity_builtin_elsets_orifield (pSim, pTess, pNodes, Mesh, + entity, res, &SimRes); + } + else nes_pproc_entity_builtin_elsets_gen (pSim, pNodes, Mesh, entity, res, members, memberqty, &SimRes); diff --git a/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin3.cpp b/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin3.cpp index d3f6d600..4cbb6a35 100644 --- a/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin3.cpp +++ b/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin3.cpp @@ -1143,3 +1143,124 @@ nes_pproc_entity_builtin_elsets_gen (struct SIM *pSim, return; } + +void +nes_pproc_entity_builtin_elsets_orifield (struct SIM *pSim, struct TESS *pTess, + struct NODES *pNodes, struct MESH *Mesh, + char *entity, char *res, struct SIMRES *pSimRes) +{ + int i, step, size, status; + double **elsetdata = ut_alloc_2d (Mesh[(*pTess).Dim].ElsetQty + 1, 4); + double *elsetdata1d = ut_alloc_1d (Mesh[(*pTess).Dim].ElsetQty + 1); + double ***evect = ut_alloc_3d (Mesh[(*pTess).Dim].ElsetQty + 1, 3, 3); + double **eval = ut_alloc_2d (Mesh[(*pTess).Dim].ElsetQty + 1, 3); + char *prev = ut_alloc_1d_char (1000); + char *filename = NULL; + FILE *file = NULL; + struct OL_SET OSet; + struct ODF Odf; + char *fct = NULL, **vars = NULL, **vals = NULL; + char *weight = NULL; + char *var = NULL; + char *thetastring = ut_alloc_1d_char (100); + struct SIMRES SimRes2; + + neut_simres_set_zero (&SimRes2); + + ut_string_function (res, &fct, &vars, &vals, &size); + + neut_sim_orispace (*pSim, &Odf, (char*) "R"); + + ut_print_progress (stdout, 0, (*pSim).StepQty + 1, (char*) "%3.0f%%", prev); + + for (step = 0; step <= (*pSim).StepQty; step++) + { + neut_sim_setstep (pSim, step); + neut_sim_simres (*pSim, (char*) "elsets", (char*) "ori", &SimRes2); + + neut_simres_setstep (pSimRes, step); + neut_simres_setstep (&SimRes2, step); + + status = neut_mesh_elsetori (Mesh[(*pTess).Dim], elsetdata); + if (status) + ut_print_message (2, 0, (char*) "Failed to read elset oris.\n"); + + file = ut_file_open ((*pSimRes).file, (char*) "W"); + + neut_mesh_elsets_olset (*pNodes, Mesh[(*pTess).Dim], elsetdata, + NULL, Mesh[(*pTess).Dim].ElsetQty, &OSet); + ut_string_string (Mesh[(*pTess).Dim].ElsetCrySym, &(OSet.crysym)); + + int user_sigma = 0; + ut_string_string ("1", &weight); + for (i = 0; i < size; i++) + { + if (!strcmp (vars[i], "theta")) + { + Odf.sigma = atof (vals[i]); + Odf.sigma *= M_PI / 180; + user_sigma = 1; + } + else if (!strcmp (vars[i], "weight")) + ut_string_string (vals[i], &weight); + else if (!strcmp (vars[i], "var")) + ut_string_string (vals[i], &var); + } + if (!user_sigma) + neut_odf_setsigma (&Odf, (char*) "avthetaeq", OSet.size, OSet.crysym); + + sprintf (thetastring, " (theta = %9.6f°) ", Odf.sigma * 180 / M_PI); + ut_print_clearline (stdout, strlen (thetastring) - 1); + printf ("%s", thetastring); + + ut_print_progress (stdout, 0, (*pSim).StepQty + 1, (char*) "%3.0f%%", prev); + + // Setting weights + for (i = 0; i < size; i++) + if (!strcmp (vars[i], "weight")) + neut_mesh_entity_expr_val (*pNodes, Mesh, pTess, + NULL, NULL, NULL, NULL, (char*) "elset", weight, + OSet.weight, NULL); + + neut_sim_simres (*pSim, (char*) "elset", var, &SimRes2); + ut_array_1d_fnscanf (SimRes2.file, elsetdata1d + 1, Mesh[(*pTess).Dim].ElsetQty, (char *) "R"); + + if (!strcmp (fct, "orifield")) + { + neut_odf_orifield_comp ((char *) "m", !strcmp (weight, "1") ? (char *) "5" : (char *) "all", &OSet, elsetdata1d + 1, &Odf); + + for (i = 0; i < Odf.odfqty; i++) + fprintf (file, REAL_PRINT_FORMAT "\n", Odf.odf[i]); + } + else if (!strcmp (fct, "orifieldn")) + { + neut_odf_orifield_comp ((char *) "n", !strcmp (weight, "1") ? (char *) "5" : (char *) "all", &OSet, elsetdata1d + 1, &Odf); + for (i = 0; i < Odf.odfnqty; i++) + fprintf (file, REAL_PRINT_FORMAT "\n", Odf.odfn[i]); + } + + ut_file_close (file, (*pSimRes).file, "W"); + + ut_print_progress (stdout, step + 1, (*pSim).StepQty + 1, (char *) "%3.0f%%", prev); + } + + neut_sim_setstep (pSim, 0); + neut_sim_addres (pSim, entity, res, NULL); + neut_sim_fprintf ((*pSim).simdir, *pSim, (char *) "W"); + + // neut_odf_free (&Odf); + ut_free_1d_char (&filename); + ut_free_3d (&evect, Mesh[(*pTess).Dim].ElsetQty + 1, 3); + ut_free_2d (&eval, Mesh[(*pTess).Dim].ElsetQty + 1); + ut_free_2d (&elsetdata, Mesh[(*pTess).Dim].ElsetQty + 1); + ut_free_1d (&elsetdata1d); + ut_free_1d_char (&prev); + ut_free_1d_char (&weight); + ut_free_1d_char (&var); + ut_free_1d_char (&thetastring); + ut_free_2d_char (&vars, size); + ut_free_2d_char (&vals, size); + neut_simres_free (&SimRes2); + + return; +} diff --git a/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin_.h b/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin_.h index 47c1a3dd..5c9bac85 100644 --- a/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin_.h +++ b/src/neper_s/nes_pproc/nes_pproc_entity/nes_pproc_entity_builtin/nes_pproc_entity_builtin_.h @@ -52,6 +52,10 @@ extern void nes_pproc_entity_builtin_elsets_odf (struct SIM *pSim, struct TESS * struct NODES *pNodes, struct MESH *Mesh, char *entity, char *res, struct SIMRES *pSimRes); +extern void nes_pproc_entity_builtin_elsets_orifield (struct SIM *pSim, struct TESS *pTess, + struct NODES *pNodes, struct MESH *Mesh, + char *entity, char *res, struct SIMRES *pSimRes); + extern void nes_pproc_entity_builtin_elsets_readodf (struct SIM *pSim, struct TESS Tess, struct MESH *Mesh, char *entity, char *res, struct SIMRES *pSimRes); diff --git a/src/neut/neut_odf/neut_odf.h b/src/neut/neut_odf/neut_odf.h index 3e2a9629..aa6737e7 100644 --- a/src/neut/neut_odf/neut_odf.h +++ b/src/neut/neut_odf/neut_odf.h @@ -24,6 +24,8 @@ extern void neut_odf_convolve (struct ODF *pOdf, char *kernel); extern void neut_odf_deconvolve (struct ODF *pOdf, char *kernel); extern void neut_odf_elt_ori (struct ODF Odf, int elt, gsl_rng *r, double *q); +extern void neut_odf_orifield_comp (char *mode, char *neigh, struct OL_SET *pOSet, + double *oridata, struct ODF *pOdf); #endif /* NEUT_ODF_H */ diff --git a/src/neut/neut_odf/neut_odf1.cpp b/src/neut/neut_odf/neut_odf1.cpp index e6e1d44f..699b6462 100644 --- a/src/neut/neut_odf/neut_odf1.cpp +++ b/src/neut/neut_odf/neut_odf1.cpp @@ -7,6 +7,10 @@ extern double neut_odf_comp_elts (char *neigh, struct OL_SET *pOSet, QCLOUD nanocloud, my_kd_tree_t *nano_index, struct ODF *pOdf, int verbosity); extern double neut_odf_comp_nodes (char *neigh, struct OL_SET *pOSet, QCLOUD nano_cloud, my_kd_tree_t *nano_index, struct ODF *pOdf, int verbosity); +extern double neut_odf_orifield_comp_elts (char *neigh, struct OL_SET *pOSet, QCLOUD nano_cloud, + my_kd_tree_t *nano_index, double *oridata, struct ODF *pOdf); +extern double neut_odf_orifield_comp_nodes (char *neigh, struct OL_SET *pOSet, QCLOUD nano_cloud, + my_kd_tree_t *nano_index, double *oridata, struct ODF *pOdf); void neut_odf_set_zero (struct ODF *pOdf) @@ -311,3 +315,24 @@ neut_odf_elt_ori (struct ODF Odf, int elt, gsl_rng *r, double *q) return; } + +void +neut_odf_orifield_comp (char *mode, char *neigh, struct OL_SET *pOSet, + double *oridata, struct ODF *pOdf) +{ + my_kd_tree_t *nano_index = nullptr; + nanoflann::SearchParams params; + QCLOUD nano_cloud; + + neut_oset_kdtree (pOSet, &nano_cloud, &nano_index); + + if (strstr (mode, "m") || strstr (mode, "n")) + neut_odf_orifield_comp_elts (neigh, pOSet, nano_cloud, nano_index, oridata, pOdf); + + if (strstr (mode, "n")) + neut_odf_orifield_comp_nodes (neigh, pOSet, nano_cloud, nano_index, oridata, pOdf); + + delete nano_index; + + return; +} diff --git a/src/neut/neut_odf/neut_odf2.cpp b/src/neut/neut_odf/neut_odf2.cpp index 8f7b749a..31a0e36d 100644 --- a/src/neut/neut_odf/neut_odf2.cpp +++ b/src/neut/neut_odf/neut_odf2.cpp @@ -11,6 +11,10 @@ extern void neut_odf_comp_elts_neigh (double *q, double cut_fact, QCLOUD qcloud, extern void neut_odf_comp_nodes_all (struct OL_SET *pOSet, double *q, struct ODF *pOdf, int id); extern void neut_odf_comp_nodes_neigh (double *q, double cut_fact, QCLOUD qcloud, my_kd_tree_t *nano_index, struct ODF *pOdf, int id); +extern void neut_odf_orifield_comp_elts_all (struct OL_SET *pOSet, double *oridata, double *q, + int id, struct ODF *pOdf); +extern void neut_odf_orifield_comp_nodes_all (struct OL_SET *pOSet, double *oridata, double *q, + int id, struct ODF *pOdf); void neut_odf_init_eltweight (struct ODF *pOdf) @@ -181,3 +185,126 @@ neut_odf_comp_nodes (char *neigh, struct OL_SET *pOSet, QCLOUD nano_cloud, return; } + +double +neut_odf_orifield_comp_elts (char *neigh, struct OL_SET *pOSet, QCLOUD nano_cloud, + my_kd_tree_t *nano_index, double *oridata, struct ODF *pOdf) +{ + int i, dim = neut_mesh_array_dim ((*pOdf).Mesh); + double *vol = ut_alloc_1d ((*pOdf).odfqty); + double avradeq; + neut_ori_n_avradeq (NULL, (*pOSet).size, (*pOSet).crysym, &avradeq); + + ut_array_1d_zero ((*pOdf).odf, (*pOdf).odfqty); + +#pragma omp parallel for schedule(dynamic) + for (i = 0; i < (*pOdf).Mesh[dim].EltQty; i++) + { + double *q = ol_q_alloc (); + double *coo = ut_alloc_1d (3); + + if (dim == 3) + neut_mesh_elt_volume ((*pOdf).Nodes, (*pOdf).Mesh[dim], i + 1, vol + i); + else + neut_mesh_elt_area ((*pOdf).Nodes, (*pOdf).Mesh[dim], i + 1, vol + i); + neut_mesh_elt_centre ((*pOdf).Nodes, (*pOdf).Mesh[dim], i + 1, coo); + + if (!strncmp ((*pOdf).gridtype, "rodrigues", 9)) + { + vol[i] *= pow (M_PI * (1 + ut_vector_scalprod (coo, coo)), -2); + ol_R_q (coo, q); + } + else if (!strncmp ((*pOdf).gridtype, "euler-bunge", 11)) + { + if (!strcmp ((*pOdf).gridunit, "degree")) + { + vol[i] *= sin (coo[1] * M_PI / 180); + ol_e_q (coo, q); + } + else + { + vol[i] *= sin (coo[1]); + ol_e_q_rad (coo, q); + } + } + else if (!strcmp ((*pOdf).gridtype, "homochoric")) + ol_homochoric_q (&((*pOdf).hfct), coo, q); + else + ut_print_message (2, 2, + "Unkown grid type `%s' (option -odfgridtype).\n", + (*pOdf).gridtype); + + if (1 || !strcmp (neigh, "all")) + neut_odf_orifield_comp_elts_all (pOSet, oridata, q, i, pOdf); + else + { + double cut_fact = atof (neigh); + neut_odf_comp_elts_neigh (q, cut_fact, nano_cloud, nano_index, pOdf, i); + } + + ol_q_free (q); + ut_free_1d (&coo); + } + + (*pOdf).odfmean = ut_array_1d_mean ((*pOdf).odf, (*pOdf).odfqty); + (*pOdf).odfsig = ut_array_1d_wstddev ((*pOdf).odf, vol, 1, (*pOdf).odfqty); + (*pOdf).odfmin = ut_array_1d_min ((*pOdf).odf, (*pOdf).odfqty); + (*pOdf).odfmax = ut_array_1d_max ((*pOdf).odf, (*pOdf).odfqty); + + ut_free_1d (&vol); + + return 1; +} + +void +neut_odf_orifield_comp_nodes (char *neigh, struct OL_SET *pOSet, QCLOUD nano_cloud, + my_kd_tree_t *nano_index, double *oridata, struct ODF *pOdf) +{ + int i; + double mean; + + ut_array_1d_zero ((*pOdf).odfn, (*pOdf).odfnqty); + + mean = 0; +#pragma omp parallel for schedule(dynamic) + for (i = 0; i < (*pOdf).odfnqty; i++) + { + double *q = ol_q_alloc (); + double *qs = ol_q_alloc (); + + if (!strncmp ((*pOdf).gridtype, "rodrigues", 9)) + ol_R_q ((*pOdf).Nodes.NodeCoo[i + 1], q); + else if (!strncmp ((*pOdf).gridtype, "euler-bunge", 11)) + { + if (!strcmp ((*pOdf).gridunit, "degree")) + ol_e_q ((*pOdf).Nodes.NodeCoo[i + 1], q); + else + ol_e_q_rad ((*pOdf).Nodes.NodeCoo[i + 1], q); + } + else if (strcmp ((*pOdf).gridtype, "homochoric")) + ol_homochoric_q (&((*pOdf).hfct), (*pOdf).Nodes.NodeCoo[i + 1], q); + else + ut_print_message (2, 2, + "Unkown grid type `%s' (option -odfgridtype).\n", + (*pOdf).gridtype); + + if (1 || !strcmp (neigh, "all")) + neut_odf_orifield_comp_nodes_all (pOSet, oridata, q, i, pOdf); + else + { + double cut_fact; + sscanf (neigh, "%lf", &cut_fact); + neut_odf_comp_nodes_neigh (q, cut_fact, nano_cloud, nano_index, pOdf, i); + } + +#pragma omp atomic + mean += (*pOdf).odfn[i]; + + ol_q_free (q); + ol_q_free (qs); + } + + mean /= (*pOdf).odfnqty; + + return; +} diff --git a/src/neut/neut_odf/neut_odf3.cpp b/src/neut/neut_odf/neut_odf3.cpp index 4030704c..9bc836b2 100644 --- a/src/neut/neut_odf/neut_odf3.cpp +++ b/src/neut/neut_odf/neut_odf3.cpp @@ -61,6 +61,92 @@ neut_odf_comp_nodes_all (struct OL_SET *pOSet, double *q, struct ODF *pOdf, int return; } +// input: q & id: orientation and id of the mesh elt +void +neut_odf_orifield_comp_elts_all (struct OL_SET *pOSet, double *oridata, double *q, + int id, struct ODF *pOdf) +{ + int i, j, n = ol_crysym_qty ((*pOSet).crysym); + double theta, tmp, *oridata_weight = ut_alloc_1d ((*pOSet).size), *qs = ol_q_alloc (); + + for (i = 0; i < (int) (*pOSet).size; i++) + { + // Computing weight (sum of gaussian kernels) + for (j = 1; j <= n; j++) + { + ol_q_crysym ((*pOSet).q[i], (*pOSet).crysym, j, qs); + ol_q_q_misori_rad (q, qs, &theta); + + tmp = 1 / sqrt (2 * M_PI * (*pOdf).sigma * (*pOdf).sigma) + * exp (-theta * theta / (2 * (*pOdf).sigma * (*pOdf).sigma)); + + oridata_weight[i] += tmp; + } + + // Multiplying by the orientation weight + oridata_weight[i] *= (*pOSet).weight[i]; + } + + // Computing weighted average + // weighting is done respective to the the previously computed weights and + // the orientation weights, if defined + + // we use (*pOdf).odf as a placeholder, as we are not actually computing an odf + (*pOdf).odf[id] = 0; + for (i = 0; i < (int) (*pOSet).size; i++) + (*pOdf).odf[id] += oridata_weight[i] * oridata[i]; + (*pOdf).odf[id] /= ut_array_1d_sum (oridata_weight, (int) (*pOSet).size); + + ol_q_free (qs); + ut_free_1d (&oridata_weight); + + return; +} + +// same as neut_odf_orifield_comp_elts_all, except that we copy to odfn[id] +// input: q & id: orientation and id of the mesh elt +void +neut_odf_orifield_comp_nodes_all (struct OL_SET *pOSet, double *oridata, double *q, + int id, struct ODF *pOdf) +{ + int i, j, n = ol_crysym_qty ((*pOSet).crysym); + double theta, tmp, *oridata_weight = ut_alloc_1d ((*pOSet).size), *qs = ol_q_alloc (); + + for (i = 0; i < (int) (*pOSet).size; i++) + { + // Computing weight (sum of gaussian kernels) + for (j = 1; j <= n; j++) + { + ol_q_crysym ((*pOSet).q[i], (*pOSet).crysym, j, qs); + ol_q_q_misori_rad (q, qs, &theta); + + tmp = 1 / sqrt (2 * M_PI * (*pOdf).sigma * (*pOdf).sigma) + * exp (-theta * theta / (2 * (*pOdf).sigma * (*pOdf).sigma)); + + oridata_weight[i] += tmp; + } + + // Multiplying by the orientation weight + oridata_weight[i] *= (*pOSet).weight[i]; + } + + // Computing weighted average + // weighting is done respective to the the previously computed weights and + // the orientation weights, if defined + + // we use (*pOdf).odfn as a placeholder, as we are not actually computing an odfn + + (*pOdf).odfn[id] = 0; + for (i = 0; i < (int) (*pOSet).size; i++) + (*pOdf).odfn[id] += oridata_weight[i] * oridata[i]; + (*pOdf).odfn[id] /= ut_array_1d_sum (oridata_weight, (int) (*pOSet).size); + + ol_q_free (qs); + ut_free_1d (&oridata_weight); + + return; +} + void neut_odf_comp_elts_neigh (double *q, double cut_fact, QCLOUD qcloud, my_kd_tree_t *nano_index, struct ODF *pOdf, int id)