Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/SnowDriftFEInit.cc

1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18#include <assert.h>
19#include <iostream>
20#include <cstdlib>
21#include <string>
22
23#include <meteoio/MeteoIO.h>
24
25#include <alpine3d/snowdrift/SnowDrift.h>
26
27using namespace mio;
28
29/**
30 * @brief Initializes system size and all data members which are required for the finite element solution
31 */
32void SnowDriftA3D::InitializeFEData( void )
33{
34nNodes = nx * ny * nz;
35nElements = (nx-1) * (ny-1) * (nz-1);
36nDOF = (nx-2) * (ny-2) * (nz-2);
37
38//LH_BC
39nDOF=nNodes;
40
41nNZ = numberOfNonzeros();
42
43nz_face.resize(6,2);
44ny_face.resize(6,2);
45nx_face.resize(6,2);
46
47nz_bar.resize(12,2);
48ny_bar.resize(12,2);
49nx_bar.resize(12,2);
50
51nz_interior.resize(1,2);
52ny_interior.resize(1,2);
53nx_interior.resize(1,2);
54
55n_corner.resize(8);
56
57qPoint.resize(3,8);
58
59classifySubdomain();
60setQuadraturePoints();
61
62//set some parameter
63//parameter to play with
64qualla = 1.0; // paramater to vary the deltak coefficient of
65// the SUPG approach
66
67theta = 1.0; // here you can vary the heigth above the point used in
68// the diffusional part of the computeDepositionFlux function
69
70sA.resize( nNZ );
71sB.resize( nNZ );
72colA.resize( nNZ );
73rowA.resize( nDOF + 1 );
74Psi.resize( nDOF );
75c.resize( nDOF );
76q.resize( nDOF );
77T.resize( nDOF );
78T00.resize( nDOF );
79
80rhs.resize( nDOF );
81precond.resize( nDOF );
82f.resize( nDOF );
83c00.resize( nNodes );
84
85q00.resize( nNodes );
86nnzA.resize( nDOF );
87adjA.resize( nNZ );
88
89//LH_BC
90gDirichlet.resize( nDOF );
91gNeumann.resize( nDOF );
92gamma.resize( nDOF );
93
94}
95
96
97/**
98 * @brief Description...
99 * this functions sets summation limits for the loops
100 * over boundary elements (faces,bars,corners), presently the assemble
101 * routines depend on these limits => fragile
102 * the values of the limits is also affected by the
103 * enumeration of the elements presently the the first
104 * element has the number 1 (instead of 0) => even more
105 * fragile
106
107 * This setting will become important for the domain
108 * decomposition, since then an element on the domain
109 * boundary is not necessarily an element of the
110 * physical boundary, and these array account for that.
111 * Note: these summation limits are related to the
112 * number of nodes, the number of elements, the number
113 * of dofs and accordingly the number of nonzero matrix
114 * elemens in the system matrix. It would be better to
115 * have one class which holds all these information
116 * consistently.
117 */
118void SnowDriftA3D::classifySubdomain( void )
119{
120// determine the summation limits for layer, depth and width
121nz_face(0,0)=1; nz_face(0,1)=(nz-3);
122nz_face(1,0)=1; nz_face(1,1)=(nz-3);
123nz_face(2,0)=1; nz_face(2,1)=(nz-3);
124nz_face(3,0)=1; nz_face(3,1)=(nz-3);
125nz_face(4,0)=0; nz_face(4,1)=0;
126nz_face(5,0)=(nz-2); nz_face(5,1)=(nz-2);
127
128ny_face(0,0)=0; ny_face(0,1)=0;
129ny_face(1,0)=1; ny_face(1,1)=(ny-3);
130ny_face(2,0)=(ny-2); ny_face(2,1)=(ny-2);
131ny_face(3,0)=1; ny_face(3,1)=(ny-3);
132ny_face(4,0)=1; ny_face(4,1)=(ny-3);
133ny_face(5,0)=1; ny_face(5,1)=(ny-3);
134
135nx_face(0,0)=1; nx_face(0,1)=(nx-3);
136nx_face(1,0)=(nx-2); nx_face(1,1)=(nx-2);
137nx_face(2,0)=1; nx_face(2,1)=(nx-3);
138nx_face(3,0)=0; nx_face(3,1)=0;
139nx_face(4,0)=1; nx_face(4,1)=(nx-3);
140nx_face(5,0)=1; nx_face(5,1)=(nx-3);
141
142
143nz_bar(0,0)=0; nz_bar(0,1)=0;
144nz_bar(1,0)=0; nz_bar(1,1)=0;
145nz_bar(2,0)=0; nz_bar(2,1)=0;
146nz_bar(3,0)=0; nz_bar(3,1)=0;
147nz_bar(4,0)=(nz-2); nz_bar(4,1)=(nz-2);
148nz_bar(5,0)=(nz-2); nz_bar(5,1)=(nz-2);
149nz_bar(6,0)=(nz-2); nz_bar(6,1)=(nz-2);
150nz_bar(7,0)=(nz-2); nz_bar(7,1)=(nz-2);
151nz_bar(8,0)=1; nz_bar(8,1)=(nz-3);
152nz_bar(9,0)=1; nz_bar(9,1)=(nz-3);
153nz_bar(10,0)=1; nz_bar(10,1)=(nz-3);
154nz_bar(11,0)=1; nz_bar(11,1)=(nz-3);
155
156ny_bar(0,0)=0; ny_bar(0,1)=0;
157ny_bar(1,0)=1; ny_bar(1,1)=(ny-3);
158ny_bar(2,0)=(ny-2); ny_bar(2,1)=(ny-2);
159ny_bar(3,0)=1; ny_bar(3,1)=(ny-3);
160ny_bar(4,0)=0; ny_bar(4,1)=0;
161ny_bar(5,0)=1; ny_bar(5,1)=(ny-3);
162ny_bar(6,0)=(ny-2); ny_bar(6,1)=(ny-2);
163ny_bar(7,0)=1; ny_bar(7,1)=(ny-3);
164ny_bar(8,0)=0; ny_bar(8,1)=0;
165ny_bar(9,0)=0 ; ny_bar(9,1)=0;
166ny_bar(10,0)=(ny-2); ny_bar(10,1)=(ny-2);
167ny_bar(11,0)=(ny-2); ny_bar(11,1)=(ny-2);
168
169nx_bar(0,0)=1; nx_bar(0,1)=(nx-3);
170nx_bar(1,0)=(nx-2); nx_bar(1,1)=(nx-2);
171nx_bar(2,0)=1; nx_bar(2,1)=(nx-3);
172nx_bar(3,0)=0; nx_bar(3,1)=0;
173nx_bar(4,0)=1; nx_bar(4,1)=(nx-3);
174nx_bar(5,0)=(nx-2); nx_bar(5,1)=(nx-2);
175nx_bar(6,0)=1; nx_bar(6,1)=(nx-3);
176nx_bar(7,0)=0; nx_bar(7,1)=0;
177nx_bar(8,0)=0; nx_bar(8,1)=0;
178nx_bar(9,0)=(nx-2); nx_bar(9,1)=(nx-2);
179nx_bar(10,0)=(nx-2); nx_bar(10,1)=(nx-2);
180nx_bar(11,0)=0; nx_bar(11,1)=0;
181
182
183n_corner[0] = 0*(nx-1)*(ny-1) + 0*(nx-1) + 0;
184n_corner[1] = 0*(nx-1)*(ny-1) + 0*(nx-1) + (nx-2);
185n_corner[2] = 0*(nx-1)*(ny-1) + (ny-2)*(nx-1) + (nx-2);
186n_corner[3] = 0*(nx-1)*(ny-1) + (ny-2)*(nx-1) + 0;
187n_corner[4] =(nz-2)*(nx-1)*(ny-1) + 0*(nx-1) + 0;
188n_corner[5] =(nz-2)*(nx-1)*(ny-1) + 0*(nx-1) + (nx-2);
189n_corner[6] =(nz-2)*(nx-1)*(ny-1) + (ny-2)*(nx-1) + (nx-2);
190n_corner[7] = (nz-2)*(nx-1)*(ny-1) + (ny-2)*(nx-1) + 0;
191
192nz_interior(0,0) = 1; nz_interior(0,1) = nz-3;
193ny_interior(0,0) = 1; ny_interior(0,1) = ny-3;
194nx_interior(0,0) = 1; nx_interior(0,1) = nx-3;
195
196}
197
198/**
199 * @brief Number of nonzero elements of the finite element system matrix
200 * Description : computes the number of nonzero elements of the finite
201 * element system matrix, assumes a rectangular domain with simple
202 * cubic connectivity of size nx ny nz and linear elements. Then the
203 * support of the basis function at each of the (nx-2)*(ny-2)*(nz-2)
204 * degrees of freedom overlaps with the support of the basis functions
205 * of a certain number of nodes. The number of these nodes depends on the
206 * location of the degree of freedom. This defines four classes of
207 * dofs: interior (= adjacent to 0 boundary nodes) face (= adjacent to
208 * one boundary node) bar (= adjacent to two boundary nodes) and
209 * corner (= adjacent to three boundary nodes) where here adjacency refers
210 * to the simple cubic lattice graph.
211 *
212 * accesses the SnowDrift variable nNZ
213 * Comments : For parallelization issues this must be in accordance
214 * with the summation limits set in SnowDriftA3D::classifySubdomain()
215 */
216int SnowDriftA3D::numberOfNonzeros()
217{
218int nnz = 27 * (nx-2) * (ny-2) * (nz-2) +
21918 * 2 * ( (nx-2) * (ny-2) + (nz-2) * (ny-2) + (nz-2) * (nx-2) ) +
22012 * 4 * ( (nx-2) + (ny-2) + (nz-2) ) +
2218 * 8;
222
223return nnz;
224}
225
226/**
227 * @brief prepare sparse matrix
228 * Description : generates the auxiliary arrays (i.e. the arrays
229 * row-pointer and column-index) for a sparse matrix of CSR
230 * format. The sparsity pattern is assumed to be given by a finite
231 * element problem on a nodal mesh which is topologically equivalent
232 * to a simple cubic lattice with N=nx*ny*nz number of nodes. For
233 * these hexahedral elements each node belongs to 8 elements. Hence,
234 * the support of the basis function of that given node overlaps with
235 * the support of 27 nodes' basis functions (including its own) which
236 * implies 27 nonzero entries per node in the interior of the
237 * domain. On the surface one has 18, 12 and 8 nonzero entries
238 * depending on the location (interior, edge, corner) From that the
239 * total number of nonzero elements can be computed. Usually, the
240 * dimensions of colA and rowA are nNZ and nDOF + 1,respectively where
241 * the last entry of rowA holds the number of nonzeros. Note: Here the
242 * sizes are nNZ + 1 for colA nDOF + 2 for rowA. This stems from the
243 * fact that originally Numerical recipes routines are used for the
244 * sparse matrix arrays in numerical recipes start with index one. In
245 * order to adopt them an the size was simply enlarged by one and the
246 * zero-index entry is unused. => fragile
247*/
248void SnowDriftA3D::prepareSparseMatrix( CIntArray& colA_loc,
249 CIntArray& rowA_loc,
250 CDoubleArray& adjA_loc)
251{
252
253// init aux variables for the loop
254int nnz = 0;
255
256rowA_loc[0]=0;
257rowA_loc[nDOF]=colA_loc.getNx();
258
259//loop over all degrees of freedom to compute the number of non-zero
260//matrix elements per dof
261for ( unsigned int iz = 0; iz < nz; iz++) {
262for ( unsigned int iy = 0; iy < ny; iy++) {
263for ( unsigned int ix = 0; ix < nx; ix++) {
264//node takes values from 1...nDOF
265const unsigned int node = iz * nx*ny + iy * nx + ix;
266nnz = 0;
267
268unsigned int diagonal = rowA_loc[node] + nnz; //HACK: to initialize with something...
269
270// loop over adjacent nodes ( adjacency defined by graph of the
271// sparsity pattern
272for ( int kz = (signed)iz-1; kz <= (signed)iz+1; kz++) {
273for ( int ky = (signed)iy-1; ky <= (signed)iy+1; ky++) {
274for ( int kx = (signed)ix-1; kx <= (signed)ix+1; kx++) {
275const unsigned int neighborNode = kz * (nx)*(ny) + ky * (nx) + kx;
276
277// if adjacent node is not a degree of freedom
278if ( ( kx < 0 || kx >= (signed)nx ) || ( ky < 0 || ky >= (signed)ny )
279|| ( kz < 0 || kz >= (signed)nz ) ) {
280//skip the node
281} else {
282//insert column index
283colA_loc[ rowA_loc[node] + nnz ] = neighborNode;
284
285//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
286//LH_DEBUG: Algorithm test BEGIN
287adjA_loc[ rowA_loc[node] + nnz ] = -1;
288
289if ( (neighborNode == node) ) {
290diagonal = rowA_loc[node] + nnz;
291}
292//LH_DEBUG: Algorithm test END
293//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
294
295nnz++;
296}
297}
298}
299}
300rowA_loc[ node + 1 ] = rowA_loc[ node ] + nnz;
301
302//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303//LH_DEBUG: Algorithm test BEGIN
304// nnzA[ node ] = nnz;
305adjA_loc[ diagonal ] = nnz+node/(nx*ny);
306// precond[ node ] = nnz+node/(nx*ny);//nnz;
307//LH_DEBUG: Algorithm test END
308//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
309}
310}
311}
312
313}
314
315/**
316 * @brief Set Bottom Boundary Condition
317 * Set the boundary condition of bottom layer +intialize concentration in rest of domain.
318 * @param c00 initial concentration/humidity at bottom
319 */
320void SnowDriftA3D::setBC_BottomLayer(CDoubleArray& var00, const param_type param)
321{
322
323if (param==CON){
324//set BC for the bottom layer ( ext. node layer number one)
325for (unsigned int iy = 0; iy < ny; iy++) {
326for (unsigned int ix = 0; ix < nx; ix++) {
327if (saltation(ix,iy) > 0) {
328var00[ iy*nx + ix ] = S_TO_H(psum.grid2D(ix,iy))/WS0 + c_salt(ix , iy);
329} else {
330var00[ iy*nx + ix ] = S_TO_H(psum.grid2D(ix,iy))/WS0;
331}
332for (unsigned int iz = 1; iz < nz; iz++) {
333var00[ nx*ny*iz + iy*nx + ix ] = S_TO_H(psum.grid2D(ix,iy))/WS0;
334}
335}
336}
337printf("BC value (0,0): %f\n", S_TO_H(psum.grid2D(0,0))/WS0);
338}else if (param==HUM){
339values_nodes_to_elements(nodes_q_ini, var00 );
340}else if (param==TEM){
341Grid3DObject tempT;
342tempT.set(nx,ny,nz,ta.cellsize,ta.llcorner);
343for (unsigned int iz = 0; iz<nz; iz++){
344for (unsigned int iy = 0; iy < ny; iy++) {
345for (unsigned int ix = 0; ix < nx; ix++) {
346tempT.grid3D(ix,iy,iz)=(nodes_Tair_ini.grid3D(ix,iy,iz)+(Cst::gravity/Constants::specific_heat_air)*nodes_z.grid3D(ix,iy,iz));
347}
348}
349}
350values_nodes_to_elements(tempT,var00);
351}
352}
353
354/**
355* @brief Set Robin Boundary Condition
356* @param aspect Location of boundary (NSWE or top/bottom)
357* @param gamma_val (double) constant to approach Dirichlet or Neumann condition
358* @param ix, iy, iz location
359* @param var00 Initial 3D field (T, q or c) that needs to be set
360* @param param Type of field (T, q or c)
361*/
362void SnowDriftA3D::setRobinBoundaryCondition(const aspect_type aspect, const double gamma_val, const int ix, const int iy, const int iz, CDoubleArray& var00, const param_type param)
363{
364const int node=nx*ny*iz+nx*iy+ix;
365
366gamma[ node ] = gamma_val;
367if (aspect==BOTTOM && (param==HUM || param==TEM)){
368gamma[node]=0.;
369}else if (aspect==BOTTOM && param==CON){
370gamma[ node ] = gamma_val;
371}
372gDirichlet[ node ] = var00[node];
373gNeumann[ node ] = 0.;
374}
375
376/**
377 * @brief Initialize system
378 * initializes vectors for the algorithm required for each time step
379 * @param colA array determining the sparsity pattern
380 * @param rowA array determining the sparsity pattern
381 * @param sA system matrix in compressed sparse row (CSR) storage format
382 * @param sB system matrix in compressed sparse row (CSR) storage format
383 * @param rhs vector which contains the right hand side of the linear system
384 * @param f source term in the advection/diffusion equation (sublimation amount)
385 * @param Psi auxiliary vector for incorporating inhomogeneous Dirichlet boundary conditions
386 * @param var variable for which equation has to be solved (Humidity, Concentration or Temperature)
387 * @param var00 boundary condition at bottom + initial condition interior domain of variable var
388 */
389
390void SnowDriftA3D::initializeSystem( CIntArray& colA_loc, CIntArray& rowA_loc, CDoubleArray& sA_loc, CDoubleArray& sB_loc, CDoubleArray& rhs_loc, CDoubleArray& f_loc, CDoubleArray& Psi_loc,CDoubleArray& var, CDoubleArray& var00,const param_type param)
391{
392for (unsigned int i = 0; i < nDOF; i++) {
393// Initialization of source/sink terms on the interior of the domain
394
395precond[i] = 0;
396Psi_loc[i] = 0;
397colA_loc[i] = 0;
398rowA_loc[i] = 0;
399rhs_loc[i] = 0;
400sA_loc[i] = 0;
401sB_loc[i] = 0;
402var[i]=0;
403var00[i]=0;
404f_loc[i] = 0;
405}
406
407rowA[nDOF] = 0;
408for ( unsigned int i = nDOF; i < nNZ; i++) {
409sA_loc[i] = 0;
410colA_loc[i] = 0;
411adjA[i] = 0;
412}
413
414setBC_BottomLayer(var00, param);
415
416for (unsigned int iz = 0; iz < nz; iz++) {
417for (unsigned int iy = 0; iy < ny; iy++) {
418for (unsigned int ix = 0; ix < nx; ix++) {
419//east
420if (( ix == 0 )
421|| ( ix == nx-1 )
422|| ( iy == ny-1 && !(ix==0 || ix==nx-1) )
423|| ( iy == 0 && !(ix==0 || ix==nx-1))
424|| ( iz == nz-1 && !(iy==0 || iy==ny-1 || ix==0 || ix==nx-1))) {
425
426//if East, West, North South or Top
427setRobinBoundaryCondition(OTHER, 1.e5, ix,iy,iz, var00, param);
428}
429//if bottom
430if ( iz == 0 && !(iy==0 || iy==ny-1 || ix==0 || ix==nx-1) ) {
431double K[3][3];
432computeDiffusionTensor(K,ix,iy, iz);
433double K_perp = K[2][2];
434setRobinBoundaryCondition(BOTTOM, K_perp,ix,iy,iz, var00, param);
435} else {
436// gamma[ node ] = 1e5;
437// gDirichlet[ node ] = 0.001;//c00[node];
438// gNeumann[ node ] = 0;//0.5 * c_salt[ ix ][ iy ];//u_perp *c00[node];
439// gDirichlet[ node ] = ( S_TO_H(sn_MdataT.nswc/WS0)
440// *( 1 + sin( 2*M_PI * ix/(1.0*(nx-1)) ) )
441// *( 1 + sin( 2*M_PI * iy/(1.0*(ny-1)) ) ) );
442}
443}//iz
444}//iy
445}//ix
446}
447
448
449//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
450// Function : resetArray
451// Authors : Henning Loewe
452// Description : set Array
453//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
454void SnowDriftA3D::resetArray( CIntArray& cda )
455{
456for (unsigned int i = 0; i < cda.getNx(); i++) {
457cda[i] = 0;
458}
459}
460
461void SnowDriftA3D::resetArray( CDoubleArray& cda )
462{
463for (unsigned int i = 0; i < cda.getNx(); i++) {
464cda[i] = 0;
465}
466}
467
468/**
469 * @brief nodes_to_elements
470 * @param nodesGrid 3Dgrid of nodes
471 * @param elementsArray array with elements
472 */
473void SnowDriftA3D::values_nodes_to_elements(const Grid3DObject& nodesGrid, CDoubleArray& elementsArray )
474{
475const int nx_grid=nodesGrid.getNx();
476const int ny_grid=nodesGrid.getNy();
477const int nxy=nx_grid*ny_grid;
478const int nz_grid=nodesGrid.getNz();
479const unsigned int Nelems=elementsArray.getNx();
480
481for (int i=0; i<(signed)Nelems; i++){
482//find the nodes for this element
483int iz = (int)floor(i/nxy);
484int iy = (int)floor(i-iz*nxy) / nx_grid;
485int ix = i-iz*nxy-iy*nx_grid;
486//elements should always be surrounded by 8 nodes, but there are some extra elements
487iz = std::min(iz, nz_grid-2);
488iy = std::min(iy, ny_grid-2);
489ix = std::min(ix, nx_grid-2);
490
491double value=0;
492unsigned int count=0;
493for (int ii=ix; ii<=ix+1; ii++){
494for (int jj=iy; jj<=iy+1; jj++){
495for (int kk=iz; kk<=iz+1; kk++){
496value += nodesGrid.grid3D(ii,jj,kk);
497count++;
498}
499}
500}
501elementsArray(i) = value/count;
502}
503}
504
505/**
506 * @brief elements_to_nodes
507 * Extract a 3D-grid of nodes from an array of elements
508 * @param nodesGrid grid of nodes that should be filled
509 * @param elementsArray array with elements
510 */
511void SnowDriftA3D::values_elements_to_nodes(Grid3DObject& nodesGrid, const CDoubleArray& elementsArray )
512{
513 const int ncols=nodesGrid.getNx();
514 const int nrows=nodesGrid.getNy();
515 const int ndepths=nodesGrid.getNz();
516 int ixmin=0;
517 int ixmax=0;
518 int iymin=0;
519 int iymax=0;
520 int izmin=0;
521 int izmax=0;
522
523 for (int ii=0; ii<=ncols-1; ii++){
524for (int jj=0; jj<=nrows-1; jj++){
525 for (int kk=0; kk<=ndepths-1; kk++){
526
527if (ii==0){
528 ixmin=ii;
529 ixmax=ii;
530}else if (ii==ncols-1){
531 ixmin=ii-1;
532 ixmax=ii-1;
533}else{
534 ixmin=ii-1;
535 ixmax=ii;
536}
537
538if (jj==0){
539 iymin=jj;
540 iymax=jj;
541}else if (jj==nrows-1){
542iymin=jj-1;
543iymax=jj-1;//iymax=jj;
544}else{
545 iymin=jj-1;
546 iymax=jj;
547}
548
549if (kk==0){
550 izmin=kk;
551 izmax=kk;
552}else if (kk==ndepths-1){
553 izmin=kk-1;
554 izmax=kk-1;
555}else{
556 izmin=kk-1;
557 izmax=kk;
558}
559
560double value = 0.;
561int count = 0;
562for ( int ix=ixmin; ix<=ixmax; ix++){
563 for (int iy = iymin; iy<=iymax; iy++) {
564for (int iz = izmin; iz<=izmax; iz++) {
565 int element = iz*(ncols)*(nrows)+iy*(ncols)+ix;
566 value += elementsArray(element);
567 count+=1;
568}
569 }
570}
571value=value/count;
572nodesGrid.grid3D(ii,jj,kk)=value;
573 }
574}
575 }
576
577}

Archive Download this file

Revision: HEAD