On Fri, Oct 20, 2017 at 12:08 PM, John Peterson <jwpeter...@gmail.com>
wrote:
>
>
> On Fri, Oct 20, 2017 at 9:53 AM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Fri, Oct 20, 2017 at 11:38 AM, John Peterson <jwpeter...@gmail.com>
>> wrote:
>>
>>>
>>>
>>> On Fri, Oct 20, 2017 at 9:19 AM, David Knezevic <
>>> david.kneze...@akselos.com> wrote:
>>>
>>>> I'm using ExodusII_IO::copy_elemental_solution to copy element data to
>>>> a new System object, and this works fine in serial using code like the
>>>> below:
>>>>
>>>> ----------------------------------------------
>>>>
>>>> Mesh mesh(comm);
>>>> ExodusII_IO exo_io(mesh);
>>>>
>>>> if (comm.rank() == 0)
>>>> {
>>>> exo_io.read(path_to_mesh_file);
>>>> }
>>>> MeshCommunication().broadcast(mesh);
>>>>
>>>> mesh.prepare_for_use();
>>>>
>>>> EquationSystems es(mesh);
>>>> ExplicitSystem& materials_system =
>>>> es->add_system<ExplicitSystem> ("data");
>>>> materials_system.add_variable ("data", CONSTANT, MONOMIAL);
>>>> es->init();
>>>>
>>>> exo_io.copy_elemental_solution(es, "data", "data");
>>>>
>>>> ----------------------------------------------
>>>>
>>>>
>>>> However, in parallel I hit the error "ERROR, ExodusII file must be
>>>> opened for reading before copying an elemental solution!" at the start of
>>>> copy_elemental_solution.
>>>>
>>>> As far as I know (e.g. based on namebased_io.C) we should only call
>>>> "read" on proc 0, so I don't see how to make this work in parallel, since
>>>> in that case exo_io is only "opened for reading" on proc 0.
>>>>
>>>> Does anyone have any suggestions on this? (John, I gather that you were
>>>> the author of most of the ExodusII reader stuff, so I thought you might
>>>> know the answer?)
>>>>
>>>
>>> Hmm... I think copy_elemental_solution() is only designed to work in
>>> serial, so could we just guard the parts the parts of that function that
>>> read from the file in "if (rank==0)"?
>>>
>>> In addition, we would only error out if the file is not open for reading
>>> on proc 0...
>>>
>>
>> OK, that makes sense. I think we could just wrap the whole method (up to
>> solution->close) in "if (rank==0)", and remove the "if" around
>> system.solution->set so that all values are set on proc 0 and then
>> communicated when we close the vector. Do you agree?
>>
>
> Yeah, that inner if-test around solution->set doesn't need to be there,
> and solution->close() is the first collective I see.
>
OK, I'll test that, and make a PR.
Question: Will the "dof_index = elem->dof_number..." line work when elem is
not owned by proc 0?
David
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
