@@ -44,199 +44,167 @@ PyObject* K_GENERATOR::getVolumeMapOfMesh(PyObject* self, PyObject* args)
4444 E_Int res;
4545 res = K_ARRAY::getFromArray3 (array, varString, f, im, jm, km, cn, eltType);
4646
47+ if (res != 1 && res != 2 )
48+ {
49+ PyErr_SetString (PyExc_TypeError,
50+ " getVolumeMap: unknown type of array." );
51+ return NULL ;
52+ }
53+
54+ E_Int api = f->getApi ();
55+ E_Int npts = f->getSize ();
4756 PyObject* tpl = NULL ;
57+ FldArrayF* f2;
4858
49- if (res == 1 || res == 2 )
59+ posx = K_ARRAY::isCoordinateXPresent (varString);
60+ posy = K_ARRAY::isCoordinateYPresent (varString);
61+ posz = K_ARRAY::isCoordinateZPresent (varString);
62+ if (posx == -1 || posy == -1 || posz == -1 )
5063 {
51- posx = K_ARRAY::isCoordinateXPresent (varString);
52- posy = K_ARRAY::isCoordinateYPresent (varString);
53- posz = K_ARRAY::isCoordinateZPresent (varString);
54- if (posx == -1 || posy == -1 || posz == -1 )
55- {
56- RELEASESHAREDB (res, array, f, cn);
57- PyErr_SetString (PyExc_ValueError,
58- " getVolumeMap: can't find coordinates in array." );
59- return NULL ;
60- }
61- posx++; posy++; posz++;
62-
63- // Coordonnees
64- E_Float* xt = f->begin (posx);
65- E_Float* yt = f->begin (posy);
66- E_Float* zt = f->begin (posz);
64+ RELEASESHAREDB (res, array, f, cn);
65+ PyErr_SetString (PyExc_ValueError,
66+ " getVolumeMap: can't find coordinates in array." );
67+ return NULL ;
68+ }
69+ posx++; posy++; posz++;
6770
68- E_Int api = f->getApi ();
69- E_Int npts = f->getSize ();
71+ // Coordonnees
72+ E_Float* xt = f->begin (posx);
73+ E_Float* yt = f->begin (posy);
74+ E_Float* zt = f->begin (posz);
7075
71- if (res == 1 ) // cas structure
76+ if (res == 1 ) // cas structure
77+ {
78+ E_Int dim = 3 ;
79+ E_Int im1 = im-1 ;
80+ E_Int jm1 = jm-1 ;
81+ E_Int km1 = km-1 ;
82+ if (im1*jm1*km1 == 0 )
7283 {
73- E_Int dim = 3 ;
74- E_Int im1 = im-1 ;
75- E_Int jm1 = jm-1 ;
76- E_Int km1 = km-1 ;
77- if (im1*jm1*km1 == 0 )
84+ if ( (im1 && (jm1 || km1)) || (jm1 && km1) ) dim = 2 ;
85+ // if ((im1*jm1)||(im1*km1)||(jm1*km1)) dim = 2;
86+ else dim = 1 ;
87+ }
88+ if (im == 1 ) im1 = 1 ;
89+ if (jm == 1 ) jm1 = 1 ;
90+ if (km == 1 ) km1 = 1 ;
91+
92+ // E_Int ncells = im1*jm1*km1;
93+ E_Int ninti = im*jm1*km1;
94+ E_Int nintj = im1*jm*km1;
95+ E_Int nintk = im1*jm1*km;
96+ E_Int nint = ninti + nintj + nintk;
97+
98+ tpl = K_ARRAY::buildArray3 (1 , " vol" , im1, jm1, km1, api);
99+ K_ARRAY::getFromArray3 (tpl, f2);
100+ E_Float* volap = f2->begin (1 );
101+
102+ // calcul du volume
103+ if (dim == 1 )
104+ K_METRIC::compSurfStruct1D (im, jm, km, xt, yt, zt, volap);
105+ else if (dim == 2 )
106+ K_METRIC::compSurfStruct2D (im, jm, km, xt, yt, zt, volap);
107+ else
108+ {
109+ FldArrayF surf (nint,3 );
110+ FldArrayF snorm (nint);
111+ FldArrayF centerInt (nint, 3 );
112+ K_METRIC::compMetricStruct (
113+ im, jm, km, ninti, nintj, nintk,
114+ xt, yt, zt,
115+ volap, surf.begin (1 ), surf.begin (2 ), surf.begin (3 ),
116+ snorm.begin (),
117+ centerInt.begin (1 ), centerInt.begin (2 ), centerInt.begin (3 ));
118+ }
119+ RELEASESHAREDS (tpl, f2);
120+ RELEASESHAREDS (array, f);
121+ return tpl;
122+ }
123+ else if (res == 2 ) // cas non-structure
124+ {
125+ // Build array contenant le volume des elements au centre des elements
126+ tpl = K_ARRAY::buildArray3 (
127+ 1 , " vol" , npts, *cn, eltType, true , api, true
128+ );
129+ K_ARRAY::getFromArray3 (tpl, f2);
130+ E_Float* vol = f2->begin (1 ); // pointeur sur le tableau de volumes
131+
132+ if (strcmp (eltType, " NGON" ) == 0 ) // Elements NGON
133+ {
134+ E_Int ierr = 0 ;
135+ if (method == 0 )
136+ ierr = K_METRIC::compVolNGon (xt, yt, zt, *cn, vol);
137+ else if (method == 1 )
138+ ierr = K_METRIC::compVolNGonImad (xt, yt, zt, *cn, vol);
139+ else
78140 {
79- if ( (im1 && (jm1 || km1)) || (jm1 && km1) ) dim = 2 ;
80- // if ((im1*jm1)||(im1*km1)||(jm1*km1)) dim = 2;
81- else dim = 1 ;
141+ PyErr_SetString (PyExc_ValueError,
142+ " getVolumeMap: wrong method (should be 0 or 1)." );
143+ RELEASESHAREDS (tpl, f2);
144+ RELEASESHAREDU (array, f, cn);
145+ return NULL ;
82146 }
83- if (im == 1 ) im1 = 1 ;
84- if (jm == 1 ) jm1 = 1 ;
85- if (km == 1 ) km1 = 1 ;
86-
87- // E_Int ncells = im1*jm1*km1;
88- E_Int ninti = im*jm1*km1;
89- E_Int nintj = im1*jm*km1;
90- E_Int nintk = im1*jm1*km;
91- E_Int nint = ninti + nintj + nintk;
92-
93- tpl = K_ARRAY::buildArray3 (1 , " vol" , im1, jm1, km1, api);
94- FldArrayF* f2;
95- K_ARRAY::getFromArray3 (tpl, f2);
96- E_Float* volap = f2->begin (1 );
97-
98- // calcul du volume
99- if (dim == 1 )
100- K_METRIC::compSurfStruct1D (im, jm, km, xt, yt, zt, volap);
101- else if (dim == 2 )
102- K_METRIC::compSurfStruct2D (im, jm, km, xt, yt, zt, volap);
103- else
147+
148+ if (ierr == 1 )
104149 {
105- FldArrayF surf (nint,3 );
106- FldArrayF snorm (nint);
107- FldArrayF centerInt (nint, 3 );
108- K_METRIC::compMetricStruct (
109- im, jm, km, ninti, nintj, nintk,
110- xt, yt, zt,
111- volap, surf.begin (1 ), surf.begin (2 ), surf.begin (3 ),
112- snorm.begin (),
113- centerInt.begin (1 ), centerInt.begin (2 ), centerInt.begin (3 ));
150+ PyErr_SetString (PyExc_ValueError,
151+ " getVolumeMap: dimension of element unknown." );
152+ RELEASESHAREDS (tpl, f2);
153+ RELEASESHAREDU (array, f, cn);
154+ return NULL ;
114155 }
115- RELEASESHAREDS (tpl, f2);
116- RELEASESHAREDS (array, f);
117156 }
118- else if (res == 2 ) // cas non-structure
157+ else // ME
119158 {
120- // Build array contenant le volume des elements au centre des elements
121- tpl = K_ARRAY::buildArray3 (1 , " vol" , npts, *cn, eltType, true , api, true );
122- FldArrayF* f2; K_ARRAY::getFromArray3 (tpl, f2);
123- E_Float* vol = f2->begin (1 ); // pointeur sur le tableau de volumes
124-
125- if (strcmp (eltType, " NGON" ) == 0 ) // Elements NGON
126- {
127- E_Int ierr = 0 ;
128- if (method == 0 )
129- ierr = K_METRIC::compVolNGon (xt, yt, zt, *cn, vol);
130- else if (method == 1 )
131- ierr = K_METRIC::compVolNGonImad (xt, yt, zt, *cn, vol);
132- else
133- {
134- PyErr_SetString (PyExc_ValueError,
135- " getVolumeMap: wrong method (should be 0 or 1)." );
136- RELEASESHAREDS (tpl, f2);
137- RELEASESHAREDU (array, f, cn);
138- return NULL ;
139- }
140-
141- if (ierr == 1 )
142- {
143- PyErr_SetString (PyExc_ValueError,
144- " getVolumeMap: dimension of element unknown." );
145- RELEASESHAREDS (tpl, f2);
146- RELEASESHAREDU (array, f, cn);
147- return NULL ;
148- }
159+ E_Int ierr;
160+ E_Int ntotFacets = 0 , ntotElts = 0 ;
161+ E_Int nc = cn->getNConnect ();
162+
163+ // Number of facets per element
164+ std::vector<E_Int> nfpe;
165+ ierr = K_CONNECT::getNFPE (nfpe, eltType, false );
166+ if (ierr != 0 )
167+ {
168+ PyErr_SetString (PyExc_TypeError,
169+ " getVolumeMap: Error computing nfpe." );
170+ RELEASESHAREDS (tpl, f2);
171+ RELEASESHAREDU (array, f, cn);
172+ return NULL ;
149173 }
150- else // ME
174+
175+ // Total number of elements/facets
176+ for (E_Int ic = 0 ; ic < nc; ic++)
151177 {
152- E_Int nc = cn->getNConnect ();
153-
154- // Get dimensionality
155- E_Int dim = K_CONNECT::getDimME (eltType);
156-
157- if (dim == 1 )
158- {
159- E_Int offset = 0 ;
160- for (E_Int ic = 0 ; ic < nc; ic++)
161- {
162- K_FLD::FldArrayI& cm = *(cn->getConnect (ic)); // TODO
163- E_Int nelts = cm.getSize ();
164- E_Int ind1, ind2;
165- E_Float dx, dy, dz;
166- for (E_Int i = 0 ; i < nelts; i++)
167- {
168- ind1 = cm (i, 1 ) - 1 ;
169- ind2 = cm (i, 2 ) - 1 ;
170- dx = xt[ind2] - xt[ind1];
171- dy = yt[ind2] - yt[ind1];
172- dz = zt[ind2] - zt[ind1];
173- vol[i+offset] = std::sqrt (dx*dx + dy*dy + dz*dz);
174- }
175- offset += nelts;
176- }
177- }
178- else if (dim == 2 )
179- {
180- // Compute total number of elements
181- E_Int ntotElts = 0 ;
182- for (E_Int ic = 0 ; ic < nc; ic++)
183- {
184- K_FLD::FldArrayI& cm = *(cn->getConnect (ic));
185- E_Int nelts = cm.getSize ();
186- ntotElts += nelts;
187- }
188-
189- // Allocate memory to store facet normals and their areas for all
190- // connectivities
191- FldArrayF snx (ntotElts), sny (ntotElts), snz (ntotElts);
192-
193- // Compute surface of elements
194- K_METRIC::compSurfUnstruct (
195- *cn, eltType, xt, yt, zt,
196- snx.begin (), sny.begin (), snz.begin (), vol
197- );
198- }
199- else if (dim == 3 )
200- {
201- // Number of facets per element
202- std::vector<E_Int> nfpe;
203- E_Int ierr = K_CONNECT::getNFPE (nfpe, eltType, false );
204- if (ierr != 0 )
205- {
206- RELEASESHAREDS (tpl, f2);
207- RELEASESHAREDU (array, f, cn);
208- return NULL ;
209- }
210-
211- // Compute total number of facets
212- E_Int ntotFacets = 0 ;
213- for (E_Int ic = 0 ; ic < nc; ic++)
214- {
215- K_FLD::FldArrayI& cm = *(cn->getConnect (ic));
216- E_Int nelts = cm.getSize ();
217- ntotFacets += nfpe[ic]*nelts;
218- }
219-
220- FldArrayF snx (ntotFacets), sny (ntotFacets), snz (ntotFacets);
221- FldArrayF surf (ntotFacets);
222- K_METRIC::compMetricUnstruct (
223- *cn, eltType,
224- xt, yt, zt,
225- snx.begin (), sny.begin (), snz.begin (), surf.begin (), vol
226- );
227- }
178+ K_FLD::FldArrayI& cm = *(cn->getConnect (ic));
179+ E_Int nelts = cm.getSize ();
180+ ntotFacets += nfpe[ic]*nelts;
181+ ntotElts += nelts;
228182 }
183+
184+ // Get dimensionality
185+ E_Int dim = K_CONNECT::getDimME (eltType);
229186
230- RELEASESHAREDS (tpl, f2);
231- RELEASESHAREDU (array, f, cn);
187+ K_FLD::FldArrayF snx (ntotFacets), sny (ntotFacets), snz (ntotFacets);
188+ K_FLD::FldArrayF surf (ntotFacets);
189+
190+ if (dim == 3 )
191+ {
192+ K_METRIC::compMetricUnstruct (
193+ *cn, eltType, xt, yt, zt,
194+ snx.begin (), sny.begin (), snz.begin (), surf.begin (), vol
195+ );
196+ }
197+ else // 1D/2D
198+ {
199+ K_METRIC::compSurfUnstruct (
200+ *cn, eltType, xt, yt, zt,
201+ snx.begin (), sny.begin (), snz.begin (), vol
202+ );
203+ }
232204 }
205+
206+ RELEASESHAREDS (tpl, f2);
207+ RELEASESHAREDU (array, f, cn);
208+ return tpl;
233209 }
234- else
235- {
236- RELEASESHAREDB (res, array, f, cn);
237- PyErr_SetString (PyExc_TypeError,
238- " getVolumeMap: unknown type of array." );
239- return NULL ;
240- }
241- return tpl;
242210}
0 commit comments