No SpecialSets for dQ0 mesh

[OlympusMonds]

Sorry for the spamming!

Here's what I do:

dim = 2

ires = 96
jres = 48

xmin, xmax = -200e3, 200e3
ymin, ymax = -160e3, 15e3

elementMesh = uw.mesh.FeMesh_Cartesian(elementType = 'Q1/dQ0', 
                                       elementRes = (ires, jres), 
                                       minCoord = (xmin, ymin), 
                                       maxCoord = (xmax, ymax))
linearMesh   = elementMesh
constantMesh = elementMesh.subMesh 

print linearMesh.specialSets.keys()
print constantMesh.specialSets.keys()

Output:

['MaxI_VertexSet', 'MinI_VertexSet', 'AllWalls', 'MinJ_VertexSet', 'MaxJ_VertexSet', 'Empty']
{}

I know it's more of cell than a vertex for dQ0, but it would be useful to have. For example, I'm trying to get the average pressure along a wall:

np.average(pressureField[constantMesh.specialSets["MaxJ_VertexSet"]])

and it fails on a KeyError.

[lmoresi]

Partly because we prefer to keep these operations at a high level through the FE variable interface. It's not obvious how to map the pressure to the v mesh unless you have access to the interpolation functions

Sent from my iPhone

On 30 Sep 2015, at 4:35 PM, Luke Mondy <notifications@github.commailto:notifications@github.com> wrote:

Sorry for the spamming!

Here's what I do:

dim = 2

ires = 96
jres = 48

xmin, xmax = -200e3, 200e3
ymin, ymax = -160e3, 15e3

elementMesh = uw.mesh.FeMesh_Cartesian(elementType = 'Q1/dQ0',
elementRes = (ires, jres),
minCoord = (xmin, ymin),
maxCoord = (xmax, ymax))
linearMesh = elementMesh
constantMesh = elementMesh.subMesh

print linearMesh.specialSets.keys()
print constantMesh.specialSets.keys()

Output:

['MaxI_VertexSet', 'MinI_VertexSet', 'AllWalls', 'MinJ_VertexSet', 'MaxJ_VertexSet', 'Empty']
{}

I know it's more of cell than a vertex for dQ0, but it would be useful to have. For example, I'm trying to get the average pressure along a wall:

np.average(pressureField[constantMesh.specialSets["MaxJ_VertexSet"]])

and it fails on a KeyError.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/7.

[OlympusMonds]

Right, I see, something more like this?

np.average(pressureField.data.reshape((jres, ires))[-1,:])  # get the average pressure at the top of the domain.

Just a note - the above is slightly confusing, as numpy uses matrix row col ordering, where as the elements are x y, so some backwards indexing is required.

Cheers!

[jmansour]

Luke you can also do a volume integral, but only over the required element. Something like this perhaps

pressure_element_integral  = uw.utils.Integral( feMesh=linearMesh,   fn=pressureField*fn.branch.conditional({ fn.input[1] < (minY + elementHeight) : 1.,
                                                                                          True: 0. } ) )
result = pressure_element_integral.evaluate() / elementHeight

This will also work in parallel. We will have surface integrals up and running shortly.

[lmoresi]

One real disadvantage of the numpy interface is having to handle communication explicitly (closely followed by loss of information about the context of the array data).

The surface integral approach is the obvious way to go. The trickiness of the implementation is a warning to me that if we don’t do it in the C code, then it will be badly butchered in a thousand different python codes !

L

On 30 September 2015 at 5:38:13 pm, jmansour (notifications@github.commailto:notifications@github.com) wrote:

This will also work in parallel. We will have surface integrals up and running shortly.

---

Professor Louis Moresi
louis.moresi@unimelb.edu.auhttp://mailto:louis.moresi@unimelb.edu.au/
(w) +61 3 8344 1217tel://8344%201217
(m) +61 4 0333 1413tel://0333%201413
(us) +1 505 349 4425tel://349%204425
www.moresi.infohttp://www.moresi.info/
www.facebook.com/underworldcodehttp://www.facebook.com/underworldcode