Dear Ranjeet,
please ask your questions to the mailing list so that other people can
profit from this discussion, too.
If you implemented your stress output similar to the way as it is done
in the UNCOUPLED poroelastic test case (which does not use multidomain),
then it is clear where the error comes from.
Currently, in the multidomain framework, a "coupling context" is
precomputed prior to element-local calculations. This is done in the
assembler before starting the assembly of element-local entries in the
global system matrix. This concept was introduced to speed up the
computations and it is somewhat similar to what you do with the geometry
and the volume variables, i.e.:
fvGeometry.bind( ... )
elemVolVars.bind( ... )
For the coupling context the call is done on the coupling manager:
couplingManager.bindCouplingContext(domainId, element, assembler)
This precomputes all variables of the other domain that are necessary
for computations on the element of your current domain. So in your
element loop of the stress computation you would have to do this call
PRIOR to the element volume variables bind() call.
We are not happy with this and are currently working on getting rid of
the necessity of this call. Until then, you will have to pass the
assembler to your stress computation routine in order to realize the
bindCouplingContext() call before computing the element stresses. Sorry
about that!
Let us know if this fixes your problem!
Best wishes,
Dennis
On 29.01.19 10:25, Timo Koch wrote:
-------- Forwarded Message --------
Subject: Re: [DuMuX] how to add stresses to output file for
multidomain-el2p problem
Date: Tue, 29 Jan 2019 19:52:41 +1100
From: Ranjeet kumar <[email protected]>
To: Timo Koch <[email protected]>
Thank you for the reply.
As per gdb & valgrind, problem seems coming from the:
const auto effPress = problem.effectivePorePressure(element,
fvGeometry, elemVolVars, fluxVarsCache); and hence from:
const auto sigma = StressType::stressTensor(problem, element,
fvGeometry, elemVolVars, fluxVarsCache);
I have attached the valgrind log file. If you could give some idea to
solve this, it will be great helpful.
On Thu, Jan 24, 2019 at 7:18 PM Timo Koch
<[email protected]
<mailto:[email protected]>> wrote:
Hi Ranjeet,
please use a debugging tool like valgrind to find where your
segmentation fault comes from. Maybe then we can help.
Timo
> Am 24.01.2019 um 08:08 schrieb Ranjeet kumar
<[email protected] <mailto:[email protected]>>:
>
> Hi All,
> I want to add stresses value to the vtu file for the el2p
problem. Based on geomechanics/poroselastic/main.cc, I tried to
implementing the same for el2p problem. But I getting error
segmentation fault error.
>
> // function to evaluate the element stresses
> template< class StressType,
> class Problem,
> class GridGeometry,
> class GridVariables,
> class SolutionVector,
> class SigmaStorage>
> void assembleElementStresses(SigmaStorage& sigmaStorage,
> SigmaStorage& effSigmaStorage,
> const Problem& problem,
> const GridGeometry& fvGridGeometry,
> const GridVariables& gridVariables,
> const SolutionVector& x)
> {
> for (const auto& element : elements(fvGridGeometry.gridView()))
> {
> auto fvGeometry = localView(fvGridGeometry);
> auto elemVolVars =
localView(gridVariables.curGridVolVars());
>
> fvGeometry.bind(element);
> elemVolVars.bind(element, fvGeometry, x);
>
> // evaluate flux variables cache at cell center
> using FluxVarsCache = typename
GridVariables::GridFluxVariablesCache::FluxVariablesCache;
> FluxVarsCache fluxVarsCache;
> fluxVarsCache.update(problem, element, fvGeometry,
elemVolVars, element.geometry().center());
>
> //get lame parameters, the pressure and compute stress
tensor
> const auto sigma = StressType::stressTensor(problem,
element, fvGeometry, elemVolVars, fluxVarsCache); // this
gives error
> const auto effSigma =
StressType::effectiveStressTensor(problem, element, fvGeometry,
elemVolVars, fluxVarsCache);
>
>
> //std::vector<int> v{1,2,3};
> const auto eIdx =
fvGridGeometry.elementMapper().index(element);
>
> // pass values into storage container
> using FVGridGeometry = typename GridVariables::GridGeometry;
> for (int dir = 0; dir <
FVGridGeometry::GridView::dimension; ++dir)
> {
> sigmaStorage[dir][eIdx] = sigma[dir];
> //sigmaStorage[dir][eIdx] = v[dir]; //sigma[dir];
> effSigmaStorage[dir][eIdx] = effSigma[dir];
> //effSigmaStorage[dir][eIdx] = v[dir]; //effSigma[dir];
> }
> }
> }
>
> In main method,
>
> //containers to store the sigma
> std::array<std::vector<ForceVector>, dim> sigmaStorage;
> std::array<std::vector<ForceVector>, dim> effSigmaStorage;
> const auto numCells =
poroMechFvGridGeometry->gridView().size(0);
> std::for_each(sigmaStorage.begin(), sigmaStorage.end(),
[numCells](auto& sigma){sigma.resize(numCells);});
> std::for_each(effSigmaStorage.begin(),
effSigmaStorage.end(), [numCells](auto&
effSigma){effSigma.resize(numCells);});
>
> for (int dir = 0; dir < dim; ++dir)
> {
> poroMechVtkWriter.addField(sigmaStorage[dir], "sigma_" +
std::to_string(dir), PoroMechVtkOutputModule::FieldType::element);
> poroMechVtkWriter.addField(sigmaStorage[dir],
"effSigma_" + std::to_string(dir),
PoroMechVtkOutputModule::FieldType::element);
> }
>
> using StressType = GetPropType<PoroMechTypeTag,
Properties::StressType>;
> assembleElementStresses<StressType>(sigmaStorage, effSigmaStorage,
> *poroMechProblem, *poroMechFvGridGeometry,
*poroMechGridVariables, x[poroMechId]);
>
> Any help is appreciated.
> Thank You!,
>
> Regards,
> Ranjeet
>
>
>
> _______________________________________________
> Dumux mailing list
> [email protected]
<mailto:[email protected]>
> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
